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Abstract 


The present study is concerned with the determination of parametric functions in nonlinear 
diffusion-dominated differential equations by an inverse technique. Examples of these functions 
are thermal conductivity and heat capacity in a conduction problem and relative permeability 
and capillary pressure in flow through an unsaturated porous medium. The parameters to be 
estimated are embedded in a partial differential equation and are constrained by the physical 
laws of nature. Hence it is to be expected that the objective function setup to estimate 
the parameters would be continuous and display monotone properties. The most appropriate 
optimization method technique needed to extract the parameters is then naturally based on 
gradient search. The conjugate gradient algorithm has been employed as the search technique 
in the present study. The gradient of the objective function requires the calculation of the 
adjoint of the differential operator governing the physical process. The use of the adjoint 
function makes the estimated parametric functions to be physically meaningful since the mirror 
image of the natural law is simultaneously enforced in the calculations. 

The objective function has been constructed in the present work by comparing the dis- 
crete measured state variables with their values corresponding to assumed property functions. 
This leads to a correction formula for the properties when the objective function is minimized. 
This procedure has been cast in the form of a computational algorithm. The algorithm is 
iterative in nature and is strated with reasonable initial guess. 

Numerical results obtained in the present work show that the inverse procedure is ca- 
pable of reasonable predictions of the unknown parameters. Larger errors are seen under 
certain conditions for small time and points located close to the boundaries. Errors are also 
observed when the process approaches the steady state. The reconstruction of the parameters 
is generally found to be sensitive to random errors in the measurement data. 
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Chapter 1 


Introduction 


When a numerical model is used to simulate a practical problem, the parameters which 
appear in terms of coefficients and/or initial and boundary conditions of the governing 
equation must be specified. The parameters are not directly measurable from the physical 
point-of-view. The inverse approach provides a way wherein measurements of state vari- 
ables are used to determine the unknown parameters. (This may include the unknown 
initial or boundary conditions). In the inverse problem of function estimation, one can 
guess the values of the parameters. Then by minimizing the differences between estimated 
and experimental values of the state variables the required parameters can be estimated. 
Therefore, the inverse problem can be formulated as an optimization problem. 

Numerical solution of the optimization problem can be classified into the following 
three categories: Gauss-Newton, gradient search, and direct search methods. In order 
to obtain the Gauss-Newton direction it is necessary to calculate the sensitivity matrix 
in each iteration of the nonlinear least squares method. If a gradient search method is 
used to solve the optimization problem, it is necessary to calculate the gradient vectors 
during each iteration. Gradient search methods are designed to avoid the calculation of 
the sensitivity matrix. Therefore, in principle, they require less computer time. However, 
more iterations may be required for convergence. Direct search methods do not require 
the calculation of either the sensitivity matrix or the gradient vector. However, the rate 
of convergence of such methods is generally slow. 

Gradient search methods have been used during optimization encountered in the 
inverse problem. During the past two decades, there has been great interest in the use of 
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the conjugate gradient method in iterative minimization procedures applied to the solu- 
tion of constrained and unconstrained problems involving linear and nonlinear equations. 
Very recently, the method has been applied to the solution of the inverse problem. The 
commonly used gradient search method utilizes the direction of negative gradient of the 
functional as the search direction, while the conjugate gradient method uses a combi- 
nation of the negative gradient of the functional and the previous descent directions for 
search. Obviously, methods that calculate each new direction of search as a part of the 
descent direction at the last iteration are inherently more powerful than those in which 
the directions are assigned in advance. 


1.1 Importance of the Problem 


Thermal properties such as thermal conductivity and heat capacity play an important 
role in the transport of heat in a wide variety of materials. This can be gauged by the 
fact that a property like thermal conductivity can vary over six orders of magnitudes in 
commonly encountered materials, and upto ten orders in certain cases. Over and above 
the large variation that is seen from one material to another, thermal properties also 
depend on temperature and composition. Thus, in a practical application, they can vary 
from point to point as well as with time. 

Experimental determination of thermal properties as a function of temperature is 
a topic of great importance. Conventional methods for obtaining thermal properties 
include steady state, periodic heating, step response and pulse inputs. A majority of 
these methods utilize the constant property assumption. The data reduction process 
then reduces to a least squares procedure employing the error between the estimated and 
measured temperatures. 

In many applications, the parameters to be estimated may change with position, 
time or the dependent variables themselves. A prominent example is establishing the 
constitutive relations for multiphase flow in a porous media. For example, the parameters 
appearing in oil-water flow through a porous formation are relative permeability of oil, 
relative permeability of water, capillary pressure and porosity. Porosity can be assumed 
to be constant for a homogeneous medium. Relative permeability of oil and water are 
functions of water saturation. Water saturation is also a function of capillary pressure, 
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namely the difference between oil and and water phase pressures. Modeling of these 
relationships is extremely important for oil reservoir simulation. The inverse approach 
provides a systematic route for determination of these constitutive relationships. 


1.2 Literature survey 


A survey of the literature reveals that various aspect of the inverse technique have been 
addressed. These have been presented below as per the following sections: 

1.2.1 Uniqueness 

The inverse problem is often ill-posed. The ill-posedness is characterized by the nonunique- 
ness and extensive dependence on data of the identified parameters. The latter stems from 
the fact that small experimental scatter will cause serious errors in the identified param- 
eters. 


Chavent. [1974] studied the uniqueness problem in connection with parameter iden- 
tification in distributed parameter systems. In the case of non-uniqueness, the identified 
parameters differ according to the initial estimate of the parameters, and there is no rea- 
son for the estimated parameters to be close to the “true” values. As a consequence, 
the responses of the model and the system will differ for inputs other than those that 
have been used for identification. Chavent studied the uniqueness problem for two situ- 
ations: (1) the case of constant parameters and (2) the case of distributed parameters in 
space. In case 1, i.e., constant parameters, there are generally more measurements than 
unknowns. This forces the inverse problem to be unique in the sense of least squares. In 
case 2, i.e., distributed parameters, if only point measurements are available, the inverse 
problem is always nonunique. The term, point measurements refers to the situation where 
measurements are made only at a limited number of locations in the spatial domain. 

The uniqueness problem in parameter estimation is intimately related to identifi- 
ability. The notion of identifiability addresses the question of whether it is all possible 
to obtain unique solutions of the inverse problem for unknown parameters of interest in 
a mathematical model, from data collected in the spatial and time domains. Kitamura 
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and Nakayiri [1977] formulated the parameter identification problem as the one-to-one 
property of the inverse problem, i.e., the one-to-one property of mapping from the space 
of system outputs to the space of parameters. However, the uniqueness of such a map- 
ping is extremely difficult to establish and often nonexistent. The authors defined the 
identifiability as follows: “We shall call an unknown parameter “identifiable” if it can be 
determined uniquely in all points of its domain by using the input-output relation of the 
system and the input-output data.” Kitamura and Nakagiri also obtained some results 
for parameter identifiability and nonidentifiability for a system characterized by linear, 
one-dimensional parabolic partial differential equation. 

An independent definition of identifiability given by Chavent [19796] is suited to the 
identification process using the output least square; error criterion. The criterion used for 
solving the inverse problem of parameter identification is said to be output-least-square- 
identifiable if and only if a unique solution of the optimization problem exists and the 
solution depends continuously on the observations. Identifiability is usually not achievable 
in the case of point measurements where data is only available at a limited number of 
locations in the spatial domain. 


1.2.2 Classification of Parameter Identification Methods 


Various techniques have been developed to solve the inverse problem of parameter identifi- 
cation. Various methods of solving inverse heat conduction problems have been discussed 
by Deck, Blackwell and Clair [1985]. Neuman [1973] classified the techniques into either 
“direct” or “indirect.” The “direct approach” treats the model parameters as dependent 
variables in a formal inverse boundary value problem. The “indirect approach” is based 
upon an output error criterion where an existing estimate of the parameters is iteratively 
improved until the model response is sufficiently close to that of the measured output. In a 
survey paper by Kubrusly [1977] on distributed parameter systems identification, the iden- 
tification procedures have been classified into three categories: (1) direct method, which 
uses optimization techniques directly to the distributed (infinite dimensional) model; (2) 
reduction to a lumped parameter system, which reduces the distributed parameter sys- 
tem to a continuous or discrete-time lumped parameter system that is described by an 
ordinary differential equation or a difference equation; and (3) reduction to an algebraic 
equation, which reduces the partial differential equation to an algebraic equation. 



1.2 Literature survey 


5 


There are two types of error criteria that have been used in the past in the formula- 
tion of the inverse problem for a distributed parameter system. Chavent ; [19796] classified 
the identification procedures into two distinctive categories based upon the error criterion 
used in the formulations. His classification is intrinsically consistent with the work of Neu- 
mara[1973]. Hence, the inverse solution methods can be classified into the following two 
categories based upon the error criterion used in the formulation of the inverse problem. 


Equation Error Criterion (Direct Method as Classified by Neuman) 

If variations of states and the derivatives (usually estimated) of those state variables are 
known over the entire domain and if the measurement and model errors arc negligible, the 
original governing equations become linear in terms of the unknown parameters. With 
the aid of boundary conditions, a direct solution for the unknown parameters may be 
possible. 

In practice, only a limited number of observations of the state variables are available. 
To formulate the inverse problem by the equation error criterion, missing data (observa- 
tions) have to be estimated by interpolation. The interpolated data in turn contain errors. 
If the interpolated data along with the observations, are substituted into the governing 
equations, some error terms will result. Such errors are called the equation errors. The 
errors are then minimized over the proper space of the parameters. It should be noted 
that approximating state variations in the entire domain using an interpolation scheme, 
without considering the statistical properties of sampling, would cause errors in the results 
of parameter identification. 

Among the available techniques we may mention the energy dissipation method 
[Nelson, 1968]; linear programming [Kleinecke, 1971]; the use of a flatness criterion [Em- 
sellem and de Marsily,, 1971]; the multiple objective decision process [ Neuman , 1973]; 
the Galerkin method [ Frind and Pinder, 1973]; the algebraic approach [Sagar et al, 
1975]; the inductive method [Nutbrown, 1975]; linear programming and quadratic- pro- 
gramming [ Hefez , 1975]; minimization of a quadratic objective function with penalty 
function [ Navarro , 1977]; and the matrix inversion method [Yeh et al, 1983]. To mini- 
mize the instability and nonuniqueness, regularity conditions are often required. 
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Output Error Criterion (Indirect Method as Classified by Neuman) 

The criterion used in this approach is generally the minimization of a “norm” of the dif- 
ference between observed and calculated states at specified observation points. The main 
advantage of this approach is that the formulation of the inverse problem is applicable 
to the situation where the number of observations is limited. It does not require differ- 
entiation of the measured data. Various optimization algorithms have been applied to 
perforin the minimization. In general, an algorithm starts from a set of initial estimates 
of the parameters and improves it in an iterative manner until the system model response 
is sullicicntly close to that of the observations. 

Among the published works in parameter identification we may mention the follow- 
ing: quasilinearization [Yeh and Tauxe, 1971; DiStcJano and Rath, 1973]; and maximum 
principle [ Lin and Yah, 1974], Yakowitz and Noren [1976]. Vermuri and Karplus [1969] 
formulated the inverse problem in terms of optimal control and solved it by a gradient 
procedure. Chen et al. [1974] and Chavent [1975] also treated the problem in an optimal 
control approach and solved it by both steepest descent method and the conjugate gra- 
dient method. Kalman filtering techniques have also been proposed in the literature for 
parameter identification [McLaughlin, 1975; Wilson et al, 1978]. Kitanidis and Vomvoris 
[1983] used the technique of maximum likelihood estimation. 

Mathematical programming techniques developed in the field of operations research 
have been utilized for solving the inverse problem of parameter identification in the field 
of petroleum engineering. Among the published works we may mention the following: 
gradient search procedures [Jacquard and Jam, 1965; Thomas et al, 1972]; decomposition 
and multilevel optimization [Haimes et al, 1968]; linear programming [Coats et al, 1970; 
Slater and Durrer, 1971; Yeh, 1975a, 6]; quadratic programming [Yeh, 1975a, b; Chang 
and Yeh, 1976]; the Gauss-Newton method [Jahns, 1966; McLaughlin, 1975]; the modified 
Gauss-Newton method [Yoon and Yeh, 1976; Yeh and Yoon, 1976; Cooley, 1977, 1982]; 
the Newton-Raphson method [Neuman and Yakowitz, 1979]; and the conjugate gradient 
method [Neuman, 1980]. 
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1.2.3 Survey of Inverse Techniques 


The heat, conduction equation is given below to illustrate typical techniques that have 
been used to solve the inverse problem. Consider the following heat conduction equation: 


d 

Ox 







subject to the initial and boundary conditions: 


(1.1) 


T(x,y,0) 

= T 0 (x,y) x,yeT = Ti+T 2 

(1.2) 

T(x, y, t) 

= 1 j (x, y, t) x, ye r. 

(1.3) 

. .OT 

A On 

= T 2 (x,y,t) x,yeT 2 

(1.4) 


For illustrational purposes, let us assume that the source term, Q is known. The parame- 
ters chosen for identification are conductivity, I( and heat capacity, C, which are assumed 
to be functions of x and t. In general, a numerical scheme is required 'to obtain solutions 
of (1.1) subject to conditions (1.2-1. 4), provided that the values of the parameters, K and 
C, are properly prescribed. Various finite-difference or finite-element methods have been 
developed for numerical simulation. In solving the inverse problem, it is essential to have 
an efficient forward solution scheme. An example is Crank-Nicolson scheme: 


+ 
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(1.5) 


The above finite-difference equations can be solvedfor T/t by an alternating direction 
method [Douglas, 1962], 
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Generalized Matrix Method Based Upon the Equation Error Criterion 


Suppose the temperature observations are available at each of the grid points and these 
observations are substituted into Equation (1.5); then the Crank-Nicolson scheme can be 
written as 


+ 

+ 
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where Ay is assumed to be equal to Arc, and 


(1.6) 


n~' nJr 2 


= ^W/'+Ty) 


'Ci., = + i<U 


To account for the lack of equality, an unknown error term 2 is added to Ecjuation 
(1.6). In practice, only a limited number of field observations is available. Interpolation 
schemes, such as cubic splines [ Ynkowitz and Noren, 1976] and kriging [Yeh at al, 1983] 
have been applied in the past to obtain values of the state at every computational grid 
associated with the numerical scheme that is based upon either finite-difference or finite 
element approximations. The error term consists of interpolation errors as well as noise 
in observations. Equation (1.6) can be simplified to 


A i — /y -j- 


1 , 2 , 


,N 


(1.7) 


with Ai the coefficient matrix, a function of T; 

Pi the vector containing values of conductivity and heat capacity at all grid points; 
N the total number of time steps; and 
bi the column vector, a function of T. 

In more compact matrix form, this becomes 


AP = b + e 


(1.8) 
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when' 


b = [b-[,%,-,b T s\ T 

e = [ef, eJ] T 

and [ Y is a transpose operator. Whether the finite difference or the finite element is 
used as the forward solution method, the resulting equation error will always have the 
form of (1.8). However, a typical finite-difference method is used to demonstrate how to 
formulate the inverse problem by the equation error criterion. The vector of unknown 
parameters, P can be determined by minimizing the equation error e. 

From (1.8), the least squares error (or residual sum of squares) can be expressed by 

e T e = (A P - b) T {A P - b) (1.9) 

Minimizing the least square error, the vector containing the parameters can be estimated 
as 


P = (. A T A)~ 1 A T b (1.10) 

where P is the estimated vector of P. The solution is highly dependent on the level of 
discretization used in the numerical solution of the governing equation. Another disad- 
vantage is that solution (1.10) is generally unstable in the presence of noise. 


Gauss-Newton Minimization Based upon the Output Error Criterion 

For modeling purposes, the objective is to determine Ii(x, y, t) from a limited number of 
observations of T(x, y, t) in the domain so that a certain cost function is optimized. If the 
classical least squares error is used to represent the output error, the objective function 
to be minimized is 

min J = [At — Ay } 1 [At — Ay) (1-11) 

K(x,y,t) 

where At is the vector of estimated temperatures based upon estimated values of param- 
eter K, and Ay is the vector of observed temperatures. 
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The Gauss-Newton algorithm has proven to be an effective algorithm to perform 
minimization. The original and modified version of the algorithm have been used by 
many researchers in the past in solving the inverse problem, e.g., Jacquard and Jain [1965]; 
Jahns [I960], Thomas at. al. [1972], Gavalas et al. [1976], Yoon and Yeh [1976], and Cooley 
[1977,1982], The popularity of the algorithm stems from the fact that it does not require 
the calculation of the Hessian matrix as is required by the Newton method and the rate of 
convergence is superior when compared to the classical gradient search procedures. The 
algorithm is basically developed from unconstrained minimization. However, constraints 
such as upper and lower bounds are easily incorporated in the algorithm with minor 
modifications. The algorithm starts with a set of initial estimates of parameters and 
converges to a local optimum. If the objective function is convex, the local optimum 
would be the global optimum. Due to the presence of noise in the observations, the 
inverse problem is usually nonconvex, and hence only a local optimum can be assured in 
the minimization. 

Let K be a vector of parameters that contains [A'i, A' 2 , • • • , A'/,]. The algorithm 
generates the following parameter sequence for an unconstrained minimization problem: 

A' n+1 = A' n — (J n P n (1.12) 

with 

A n P n = g n (1.13) 


where 

A n = [J K (I< n )} r [J K {K n )],(L x A); 

9 n = [Jk{K ,1 )] t [At — Ay], (L x 1); 

Jk Jacobian matrix of temperature with respect to A', M x L; 
(J n step size, (scalar); 

P n Gauss-Newton direction vector, (L x 1); 

M number of observations; 

L parameter dimension. 
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The elements of Jacobian matrix are represented by the* sensitivity coefficients, 

\ 


Jk 


\ 


L>Tl 

OTi 

0'jx_ 

OK i 

ok i 

OK L 

or> 2 

OT'2 

O' I 2 

QIC i 

OK 2 

QK L 

§Im 

QTm 

. . . dr M- 

dKi 

ok 2 

9K l 


(1.14) 


where M is the total number of observations, and L is the total number of parameters. 


The stop size (i n are determined by a quadratic interpolation scheme such that 
J(I \ n 1 ') < ./(/v’ n ), or simply by a trial-and-error procedure. Occasionally, the direction 
matrix [Jk 1 -h<\ become ill-conditioned. As stated earlier, the original Gauss-Newton 
algorithm does not handle constraints. 


In solving the inverse problem, we need to calculate the above sensitivity matrix in 
each iteration. 


Conjugate Gradient Method 


Neuman [ 1 980) developed an efficient, conjugate gradient algorithm for performing min- 
imization of the objective function. He extended the variational method developed by 
Chavent [1975] for calculating the gradient of the functional with respect to the param- 
eter. The conjugate gradient method uses a combination of negative gradient of the 
functional and previous descent direction as the latest search direction. In this case, the 
following iterative process is used for the estimation of K by minimizing the objective 
function (1.11) 

K n+ 1 = K n - (1.15) 

where (3% is step size for K in going from iteration n to iteration n+1, and is direction 
of descent (i.e., search direction) for I(. It is given by 

n = -n + (no) 


which is the conjugation of (a) the gradient direction J' n K at iteration n and (b) the 

for K . The conjugate coefficient, is given by 

_ iun) 2 dx 


direction of descent 1 


fU(-n-') 2 <ix 


with 


v 


K 


0 


(1.17) 
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In this optimization method, the solution converges rapidly and is not as sensitive to the 
measurement errors as compared to Gauss-Newton or other conventional optimization 
technique. More recently, the method has been applied to the solution of the inverse 
heat conduction problem. Huang and Ozisik [1992] applied this optimization technique 
for determination of wall heat flux in laminar flow and contact conductance during metal 
casting. Huang and Yan [1995] employed this method for determination of conductivity 
and heat capacity. 


1.3 Scope of the Present Work 


The present study is concerned with the determination of parameters that arise in nonlin- 
ear diffusion-dominated problems. The parameters are not constant values; they depend 
intricately on the dependent, variables themselves. They can be identified with physical 
properties such as thermal conductivity and heat capacity in heat conduction applica- 
tions. They can also be interpreted as constitutive relationships required to model flow 

in an unsaturated porous medium. 

1 

The parameter estimation procedure adopted in the present work follows a broad- 
based ‘inverse technique’. This technique is based on the minimization of a suitable 
cost function using a. gradient search algorithm. This search is augmented by using the 
conjugate gradient, method. Special features of the technique used are: 


1. the use of the adjoint operator of the governing partial differential equation to 
construct the derivative of the cost function and 

2. the use of sensitivity functions to identify clearly the portions of the measured data 
that seriously influence the predicted results. 



Chapter 2 


Mathematical Formulation 


The mathematical procedure of the inverse method based on gradient search is explained 
for a variety of steady and unsteady problems in the present chapter. The original algo- 
rithm has been adapted from the work C. Ii. Huang and Yan( 1995) and has been applied 
to steady and unsteady diffusion-like problems. The individual steps of the algorithm can 
be stated as follows: 

1. Assume the property functions, namely the parameters to be estimated. 

2. Solve the adjoint problem and obtain adjoint variables. 

3. Compute the gradients of the objective function. 

4. Compute the conjugate coefficients. 

5. Compute the directions of descent. 

G. Solve the sensitivity problems to obtain sensitivity functions. 

7. Compute search step sizes. 

8. Update the parameters to be estimated. 

9. Repeat the above computational procedure until the convergence criterion is satis- 
fied. 
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2.1 Definition of Direct and Inverse Problems 


The governing equation of a physical problem can be symbolically represented as: 

L{u,p; i) = 0 (2.1) 

where 


u = 

{u u u 2 • • -n nu ) r 

(2.2) 

p = 

(Pi 1 P‘2 

■Pn P f 

(2.3) 

L = 

(L u L 2 - 

. . T \T 
nu ) 

(2.4) 


a; is a spatial coordinate, t is time, {u\, u 2 • • • u nu ) are nu state variables and (pi,p 2 • • • p np ) 
are np parameters. The set (L x , L 2 , • • • , L nu ) represents linear/nonlinear differential oper- 
ators acting on the state variables (ui,u 2) - • • , u nu ). The initial and boundary conditions 
needed to solve Equation (2.1) are 


u 

= /o 

when t = t 0 

(2.5) 

Lub 

— fl 

Oil Fjjb 

(2.6) 

Lib 

= S 2 

on Tib 

(2.7) 


Luji and L/, ;J are vector operators representing boundary conditions on the surfaces IVb 
and V [jj respectively. 

2.1.1 Direct Problem 

The direct problem is concerned with integrating the operator L, i.e. the estimation of the 
unknown state u when parameter p and all initial and boundary conditions are prescribed. 


2.1.2 Inverse Problem 

The inverse problem carries out the estimation of unknown parameters p or initial or 
boundary conditions when continuous or discrete observations of state u are given. 
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2.1.3 Sensitivity Analysis and Sampling 


As discussed by Knopm,an and Voss [1987], the accuracy with which parameters can be 
estimated depends on the points in space and time at which the data for the inverse 
problem is collected. The parameter estimates will be accurate when the data for the 
inverse problem is obtained by sampling at points where the state variables have the 
highest sensitivity to the parameters. Examination of sensitivities is the starting point for 
designing the sampling experiment. Of interest is the magnitude of the sensitivity rather 
than its sign. Therefore the distribution in space and time, of the absolute sensitivities 


Si = 


du(x , t]p) 

0]>i 


( 2 . 8 ) 


is investigated. Sensitivity coefficients Si are to be computed for all the parameters pi 


2.2 Steady State Heat Conduction 

2.2.1 Direct Problem 


The following steady-state problem is first considered. A slab of thickness L is kept at two 
fixed temperatures at its two boundaries. The surface at x = 0 is kept at a temperature 
7} while the other boundary at x = L is kept at T r . The governing equation for steady 
heat conduction in the slab is given by 


d_ 

dx 


*< f >S 


= o 


T = To at x = 0 

T = 7} at x = L 


(2.9) 

( 2 . 10 ) 

( 2 . 11 ) 


With the following dimensionless quantities 


G © © © © © 

x=o 


© 

X=1 


Figure 2.1: Thermocouple arrangement for measurements at M points 
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we get the following dimensionless form 


d 

dx 


K(T) 


dT 

dx 


0 


T - T„ 
T = 7’, 


in 


for 

for 


0 < < 1 

(2.12) 

. 7 : = 0 

(2.13) 

X = 1 

(2.14) 


The direct problem is concerned with the determination of temperature in the slab when 
thermal conductivity and the boundary conditions at x — 0 and x = 1 are known. 


2.2.2 Inverse Problem 

For the inverse problem, thermal conductivity K(T) is regarded as the unknown, but 
other quantities in Equations (2.12-2.14) are known. In addition, temperature readings 
taken at some appropriate locations are considered available. Referring to Figure 1, it 
is assumed that M sensors are used to record the temperature information to identify 
K{T) in t,h(' inverse calculations. Let the temperature readings taken by these sensors 
be denoted by ) \ (:r * ) = 1 \,i = 1, • • • M where i — 1 and M always correspond to x = 0 
and 1 (i.e. boundary measurements) respectively. Then the inverse problem can be 
stated as follows: by utilizing the measured temperature data, Iq, estimate the unknown 
temperature-dependent, thermal property, K{T). 

In the inverse calculations, the measured temperatures are known either from nu- 
merical simulation or from real experiments. Once the temperature field is known, there 
exist some unknown but fixed thermal property that its value (a number) K(x), at any 
position, x , must satisfy the Fourier equation to give this known temperature distribution. 
Therefore, in the inverse problem of function estimation, one can replace K{T) by K(x). 
By using the minimization procedure of the objective function .7, the function K(x ) (and 
then K(T )) can be determined. The inverse problem is generated by requiring that the 
following functional be minimized: 

M 

J{K(T )] = ./[/I'M] = £ PH*) - 5H.T,)] 2 (2.15) 

1=1 

Here, 7) arc the estimated temperatures in the slab at the locations x = X{. These 
quantities are determined from the solution of the direct problem (2.12-2.14) using an 
estimated K(x) for the exact K(x). The measured temperatures at x = X{ are denoted 
by Y t . 
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2.2.3 Conjugate Gradient Method 


The following iterative process based on the conjugate gradient method is now used for 
the estimation of Ii(x) by minimizing the above functional J[K{x)\. 

I< n+ \x) = K n (x) - P n K P%{.x) (2.16) 


where is the step size for K in going from the iteration number n to n + 1 and P% is 
the direction of descent (i.e. search direction) for K. It is given by 

^w = r^(i)+</’r , w (2.i7) 


This is the con jugation of the gradient direction at iteration n with respect to the 
direction of descent P}\~ x at iteration number n — 1. The conjugate coefficient is given by 
the recursion formula 


'K 


lUoUr'Vdx 


with v 


K 


0 


(2.18) 


We note that when v = 0 for any n, in Equation (2.18), the direction of descent PkW 
becomes the gradient direction, i.e. the ‘steepest-descent’ method is obtained. 


To perform the iterations according to Equation (2.16), we need to compute the 
step sizes /f ; (- and the gradient of the functional J ,n K . In order to develop expressions for 
the determination of these quantities, the ‘sensitivity problem’ and ‘adjoint problem’ are 
constructed as described below. 


2.2.4 Sensitivity Analysis and Search Step Size 

In order to derive the sensitivity problem, it is necessary to perturb K. It is assumed 
that when K(x) undergoes a variation AA’(r), T(:r) is perturbed by T + Al\- Then 
replacing k by k + A k and T by T + A Tk in the direct problem, subtracting from the 
resulting expressions and neglecting the second-order terms, the following problem for the 
sensitivity function A Tk is obtained 


* + * 
ax [ ax J ax 

k^l=o 

ax 

(2.19) 

AT k = 0 

at x = 0 

(2.20) 

ATk = 0 

at. x = 1 

(2.21) 
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This system can be solved as a direct problem for A7 V . The functional (2.15) for iteration 
n + 1 then becomes 


./(£;» + >) = t[Ti{k n -p n k PZ)-Yi\ 2 

i= 1 

= E f ), - Uj 2 (2.22) 

i=l 

= E[Ti(k n ) - 0U^T£), - Yif 

In this derivation, it is assumed that A K = Pfi and A T%(xi) = (Pk^)i 

The step sizes /!(.' are calculated in such a way that ./( k n+l ) is a minimum. Thus 


leading to 


dJ{K n+l ) 

m 


= 0 


ft = ^ (2-23) 

£[(A WiY 

7—1 


The sensitivity function AT;" required in this formula is calculated from Equations 2.19- 
2.21. lb apply Equation 2.16, the gradient of the objective function J is needed. This is 
calculated as follows: 


2.2.5 Adjoint Problem and Gradient Equation 


To derive the adjoint problem for K{x), Equation (2.12) is multiplied by the adjoint 
function A(:/;) and the resulting expression is integrated over the space domain. Then the 
result is added to the right-hand side of equation (2.15) to yield the following expression 
for the functional J[K(x)}: 


l 

./[/,»]= / 


A 




M 


dx + Y.W-Ytf 


(2.24) 


The functional J[K(T)] is discrete due to the last term of Equation (2.24). To make this 
functional continuous we introduce the Dirac delta function. Then the functional can be 
written as 




1 


= /a 

'r = n 


iL 

dx 



(£jx) \ 
dx J 


1 m 

dx + £[r - 

Jx=0 i = 1 


Xi)dx 


(2.25) 
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which is continuous in the sense of theory of distributions. 


The variation of the functional A due to variation of K, is obtained by perturbing 
T by T + A '1'k and K by K 4- A K in Equation (2.25), subtracting from the resulting 
expression the original Equation (2.25) and neglecting the higher-order terms. We thus 
obtain 

' I d ( (U/\T,A \ d ( d'l'\ I 

dx 


A J K = j X 

.T - 0 


d d f dT' 


r\ m 

+2 / £[T - Y]AT k 8(x - Xi)dx 

Jx ~ n i = i 


(2.26) 


Or A J K = 


XK(x) 


d(AT K ) 

dx 


J 1=0 


1 

- / K(x) 


x~0 


(IX 

dx 


d(AT K ) 


dx 


dx 


+ 


. . rr dT 

1 1 
f 

dA dT' 

AAA — 

- / A A 


dx 

n J 

x=0 x=0 

dx dx 


dx + 2[{T - Y)AT k ] 


X ~1 


r l At— 1 

+2[(T - Y)AT k ] x=0 + 2 £[T- Y}AT k S(x - Xi )dx 

Ax=0 i=2 


Or A J K = 


A 

dx 


- 1 

dAl 

J X“() 



l l 

+ 

1=0 x=0 


/ A7 « 


d .dA N 
dx [ K ^dx. 


dx 


+ 


dT 

1 1 
f 

"dA dT 

AAA — 

- / A/{ 

d.r_ 

J 

x=0 x=() 

dx dx 


dx + 2[(T - Y) AT K ] 


X = 1 


-1 M - 1 

+2[(T - Y) AT*] I=0 + 2 J2 i T ~ Y]AT k 8(x - x { )dx (2.27) 

Jx=0 • o 


rl At— 1 

E 

i=2 

From sensitivity Equation (2.19)we see that ATk vanishes at the two boundaries. Hence 
Equation (2.27) reduces to 

v 1 1 1 

4- f AT, I 

dx 


AJ k 


AA-M™ 


dx 


+ J A T k 

x =° x=0 


d ( d\\ At — 1 

*(*)£) +Z[T-Y)5(x-x i ) 


dx 
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r d.r] 

1 1 

f . . 

d\ dT 

+ 

A A K — 

- / AA 

fix dx 


fix 

*-■» *-o 



dx 


(2.28) 


Since A .!/< depends only on A7\' but, does not, depend on (A'J')k, the integrands containing 
A Tk should be zero. Hence we have 


d ( d\\ M ~ l 

h + ^ ^ “ Xi) = 0 


d. 

dx 


A = 0 at x = 0 and x — 1 
Equation (2.28) thus reduces to 


(2.29) 

(2.30) 


1 

A J K = - I A K{x) 


T. — 0 


dAdT_ 
dx dx 


dx 


(2.31) 


From the definition of A Jk we have 


1 

A J K = I J' K {x)AI<{x)dx 


(2.32) 


.T = 0 


Comparing Hqua. turns (2.31) and (2.32) we get the following expression for the gradient 
J' K {x) of the functional J: 


j> ( x \ _ 

' l<( ’ dxdx 


(2.33) 


where T is calculated from direct problem (2.12-2.14) and A from the adjoint Equation 
(2.29-2.30). 


2.2.6 Discrepancy Principle for Stopping Criteria 

If the problem involves no measurement errors, the traditional convergence check specified 
as 


\J\ < e (2.34) 

where e is a small specified number, can be used as the stopping criterion. However, the 
observed temperature information contains measurement errors; as a result, the inverse 
solution will tend to approach the perturbed input data, and the solution will exhibit 
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oscillatory behavior as the number of iterations is increased. Computational experience 
has shown that it is useful to use the discrepancy principle for terminating the iteration 
process in the regular method. Assuming T(xi) — Y{x,i) = a, the discrepancy principle 
that establishes the value of the stopping criterion e can be obtained from Equation (2.15) 
as 

i 

J a 2 dx = ( 2 (2.35) 

o 

where a is the standard deviation of the measurement error. Then the stopping criterion 
is taken as 

|J| < r 2 (2.36) 

where e is determined from Equation (2.35). In the present analysis the following stopping 
criterion is used 

| J n+1 - J n \ < e (2.37) 

The iterative procedure for calculating the function K {x) can be summarized as follows: 

1. Assume the form of K{x)\ a constant function is a useful starting point. 

2. Generate the corresponding temperature field T(.r), by solving Equations (2.12- 
2.14). 

3. Solve the adjoint problem, Equations (2.29-2.30), and obtain adjoint variable A(.x). 

4. Compute J' k from Equation (2.33). 

5. Calculate the conjugate coefficient v K using Equation (2.18). 

6. Estimate the direction of descent P/<-, from Equation (2.17). 

7. Solve the sensitivity problem, Equations (2.19- 2.21), and obtain A T K . 

8. Compute step size fl K using Equation (2.23). 

9. Estimate K using Equation (2.16). 

10. Repeat the above calculation procedure until the discrepancy principle is satisfied. 

At convergence, J — y 0, T(x j) will be close to the measured data Ij and so the property 
K{T) can be recovered by correlating K(x) and T{x). 
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2.3 Transient Heat Conduction 


The procedure presented above for steady heat conduction is now extended to include 
transient heat conduction. 


2.3.1 Flux-Flux Boundary Conditions 

Direct Problem 


We consider the following unsteady heat conduction problem. A slab of thickness L is 
initially at temperature T 0 . The boundary at x = 0 is subjected to a constant heat flux 
qi while the other boundary at x = L is insulated. The governing differential equation of 
the physical problem can be written as 


dx 



= pClT) 

dT(x,t ) . - 

— )r^- in 0 < x < L 

ot 

(2.38) 

OX 

= ?I 

at x — 0 

(2.39) 

K( T) aT ^ 

OX 

“ Qr 

at x = L 

(2.40) 

n*j) 

- to 

for /, = 0 

(2,11) 


With the following dimensionless variables 


x 

V 


T=t, To = 5, K = 


T t 


t = 


K r t 

pC r L 2 ’ 


Qi = 


Ui 

KrV 


Qr — 


T r 
Ur 
KrT' 


c = 


K 

K r 

C_ 

C T 


we obtain dimensionless form of Equations (2.38 -2.41) as follows: 

_a_ 

dx 


■ 3r(M)\ 

OX J 

= C(T ) 

8T(x,t) . 

— — m 0 < x < 1 
ot 

(2.42) 

K(T) dT M 

ox 

= ® 

at x — 0 

(2.43) 

_ A sr^.i) 

OX 

= Qr 

at x = 1 

(2.44) 

T{x,t) 

= T 0 

for t = 0 

(2.45) 
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The overbar and the subscript V denote dimensional and reference quantities, respec- 
tively. 

There are two unknown parameters, K[T) and C(T) to be determined from a single 
differential equation. Hence there is a need to calculate sensitivity functions, adjoint 
functions and gradient equations for both K and C. 


Sensitivity Problem and Step Size 


The sensitivity function A Tj< is defined as the change in estimated temperatures at dif- 
ferent positions and times due to a change A K in K(x,t). From (2.42-2.45) we get the 
equations governing the sensitivity function AT*- as follows: 


d_ 

Ox 




dA T K (x,t) \ 
dx ) 


+ 


- K{x,t ) 


dA T K (x,t) 

dx 

d&T K {x,t) 

dx 

ATk(x, t) 


5ATkOM) 
dx 

A K (x, t ) — — for x — 0 

dx 


dT(x,t) \ 
dx ) 

for 0 < x < 1 


A A” (x, t) for x = 1 

dx 

0 for t = 0 


(2.46) 

(2.47) 

(2.48) 

(2.49) 


Similarly, the sensitivity function A T c corresponding to the property C is given by 




dx 


dx 


C (glS + A C(x, t) for 0 < a; < 1 (2.50) 


dx 


dx 


dAT c {x,t) 

1 =0 for x = 0 

dx 

= o torI = 1 

dx 

A T c (x, t) = 0 for t = 0 


(2.51) 

(2.52) 

(2.53) 


The solution of the inverse problem i.e.A' (.x, t) andC(x, t) is to be obtained in such a way 
that the difference between estimated results and experimental results is minimized. Here 
we want to minimize the sum of squares of differences between the estimated results and 
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experimental l( ’ sl, ^ ,s **■*' dillerent points of the slal) over a time period if. Therefore, the 
objective functional can be written as 


J[K(T),C(T)\ 


1 m 

/ EWJt.c) 

t - 0 1=1 


Yi(xi,t)] 2 dt 


(2.54) 


where, Yi(xi,~t) * s ^ ie temperature reading taken with a sensor placed at x = x* at time 
t M sensors are used to record temperature history over a time tf. All the sensors are 
assumed be uniformly placed. 

The following iterative process based on the conjugate gradient method is now proposed 
to update the properties: 

A n+1 (x, /.) = f( n (x, t) - (1%PZ(x, t.) (2.55) 

C n+l (x,t) = C n (x,t) - f%P£{x,t) (2.56) 

The directions of descent for I< and C are as follows 


P£(x,t) = J ,n K (x,t) + v n K P%-\x,t) (2.57) 

PcM = J' n c(x,t) + vlPZ~\x,t) (2.58) 


The sensitivity coefficients for K and C are given by 


v'l = 


/ f (J'l-) 2 dxdt 

t= 0 x=0 

K ~ TT~i 

5 f (JT') 2 dxdt 

t = 0 31=0 

I I (d’lfdxdt 

t=0 x=(J 

f f ( J'c~ l ) 2 dxdt 


t= 0 1=0 

The functional J(I( n+ \ G m+1 ) for iteration n + 1 is obtained as follows 

t j M 

J(K n+l ,C n+l ) = J Y,i T i(K n - p n K P", c n - PZP2) - Yi\ 2 


dt 


(2.59) 


(2. GO) 


(2.61) 


*=0 


4 = 1 


Expanding the above expression of the functional using Taylor series and neglecting 
higher-order terms we have 


7 Af 

j(I< n+ \C n+ ')= J 5Z[7f(/<' n , C n ) — /3^(AT^)j - /3c(AVc)i - Yi 

1-0 i=1 


l dt 


(2.62) 
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here the following assumptions are made. 


DT n 


AT n — p n *11 

A1 °~ Ic oc 

The step sizes /!)(. and p\) are to be estimated in such a way that the objective function 
is minimized at n + 1 iteration. Thus we have 


J(K n +\C n+l ) 

dP n K 

J(K n+] ,C n+l ) 

d(% 


= 0 

= 0 


(2.63) 

(2.64) 


which produce the following two linear equations for /?£• and ff ( ) . 


where 


OuPk + a nPc = 
“21 P'k + 0 , 22 Pc ~ ^2 


“n 

“12 

“22 

“21 

bi 


p dVl 

I 

=0 i=1 
tf r M 

4 »■ 

t f m 

/ EKATS)]?* 


*/ i 

t=0 2=1 
“12 
l l M 


A iVi 

= / £( 7 ? - U)( ATjJJj it 

l.~ (1 * =l 

'f M 

= J 'Z(T?-Y i ){AT c n )idt 


1=0 «=1 


(2.65) 

( 2 . 66 ) 

(2.67) 

( 2 . 68 ) 

(2.69) 

(2.70) 

(2.71) 

(2.72) 


From Equations (2.65-2.65) we get the following simplified expressions of step sizes /?£ 
and fig. 


nn _ (frl“22 ~ hon) 

K (“ll“22 ~ “12) 
on _ {hon — hon) 
C (“11 “22 — “12) 


(2.73) 

(2.74) 
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The gradient, of the objective function is calculated using the adjoint problem approach 
as follows. 


Adjoint Problem and Gradient Equation 


To derive the adjoint problem for I((x,t ), Equation (2.42) is multiplied by the Lagrange 
multiplier \{x , t.) and the resulting expression is integrated over the time and space do- 
mains. TIkmi the result is added to the right-hand side of Equation (2.54). Hence the 
resulting expression of the functional J[K(x,t),C(x,t)\ becomes 



dT(x,t) \ 
dx ) 


C(T) ! EhJl 


dxdt 


l 111 

+ j TiM-Yfaut)]* 


dt. 


(2.75) 


t=o *=> 


The variation A is obtained by perturbing T by A Tr in Equation (2.75), subtracting 
from the resulting expression the original Equation (2.75) and neglecting the higher-order 
terms. We thus obtain 

AJ * = / / A [s + h ( A * f ) - c ^] 

t = 0 x=0 L X \ / J 


Or A J K — 


h i 


+2 J J J2l T ~ Y\bT K &(x - Xi)dxdt 


•f ** 

t—0 x=0 


d(A T K Y 


t—0 
l f r 


XK{x, t) ~ v dt 

0X 4=0 x=0 


4 - 


/ 

!.=() 

I 

t = 0 

I [C\AT K ] , t L 0 dx+ j j A T K d -^-dxdt 

x= 0 t = o x=0 J ' 


-//**■>£ 


3(A Tk) 


dx 


dxdt 


OT 

AAA'A- 

ax 


dt — J J AK 

1=0 4=0 x=0 

'f 1 


d\dT_ 
dx dx 


(2.7C) 
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+ 2 f \(T - Y)AT K ] s=l dt + 2 J[(T- Y)AT K ] x=0 dt 

t=n i-o 

'/ )■ « > 

+ 2 J J Y.\ T - Y \^ T K8{x-x i )dxdt 


t= 0 x=0 


Or A J K = 


l 

/ 


XK(x, t) 


0(AT K ) 


1-0 
*! r 


<9x 


'/ 1 


- / 


/ o 
l 


AT K I<(x,t) 


OX 

Ox 


dt+ J j at* 

Jx=0 t=Ox=0 

<1 


d (j'< *\ dX ' 

& A(M) & 


dxdt 


dt+ / AAA" 


37=0 t=o 

'/ i 


or 


nl 1 '/ J; 

; „*-/ / 

J x=0 t=o I=0 


a# 


3A3T 
3x dx 


dxdt 


- J [CXAT K } l t L 0 dx + J J AT K ^^-dxdt + 2 j [(T - Y)AT K ] x =idt 

x—0 t = 0 x=0 t—0 

t f t f 1 M-l 

+ 2 J[(T-Y)AT K } x=0 dt + 2 J j £[T-F]A T K 8{x - xjdxdt 


t—0 x— 0 


Or A .7 k = 


'/ i 


/ / Ar * [s ^ + 2 ' Ij) 

/ ().T™.() L \ / 1-2 

+ j A T k \K( x ,t)^ + 2(T-Y) 


dxdt 


t = o 
*/ 


/ax. 


i=0 

1 


tf(*,t)f; - 2(T - Y) 


dt 


J x = l 


f \CXAT K \)L {) dx + f A 

/•* 0 7,=0 


./• 0 
'/ l 


1 r K{XJ) ^ +AK f' 1 

dx dx 


dt 


x ~ 0 


/ / A * 


3A3T 
dx dx 


dxdt 


(2.77) 


£=0 x— 0 

From the sensitivity problem at the two boundaries we get 

*(.. = 0 

a:r a:r 


(2.78) 
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Note that. A .7^ is not. a function A !)<. Therefore the integrands containing ATk s ' l0U 
be zero. We thus get 

(2.79) 


d_ 

dx 


Q(CX) 

dx dt 


M~ 1 

E 

i=2 


K(x, Oif + + 2 ^ [T - Y}6(x - Xi ) = 0 


r\ \ 

K{x,t)? r + 2{T-Y)=0 at. x = 0 

OX 


OX 

K{x, t) — 2(T — Y) = 0 at x = l 

ox 


A = 0 at, /. = /,, 


Equation (2.77) reduces to 


l f l 

A Jk = ~ f I AA' 


t= 0 i=0 


9x <9x 


dxdt 


( 2 . 80 ) 

( 2 . 81 ) 

( 2 . 82 ) 

(2.83) 


The expression for the variation of the functional A Jc due to C is derived next. FoH° w111 ® 
the same procedure as for property Ii we get 

'/ i 


A J c 


A 

f. = 0 ,T = 0 

'/ l 





AC°X-C°-^ 

dt ) dt 


dxdt 


+ 2 J I [T - Y]AT c 6{x - Xi)dxdt 


(2.84) 


i=l)i-0 i - 1 


Or A. 7 C = 




/ 


4=0 
t S 1 


AC 


<= 0 1=0 


A 


dx 

dT 

dt 


l f i 


dx 

dx 


+ 


i 

r 

cl(CA) 

J 

=0 

dt 


dt- f [ K(x,t) 

t=0 x=0 

1 

dxdt - f [CXATcftLodx 
£=0 

ATcdxdt 4- 2 j [(T-r)A7b]x= 


9(AT c ) 


dx 


dxdt 


dt 


t = o 

+ 2 J [(T-Y)AT c ) x =odt + 2 j j ^[T - Y}AT c 6(x - Xi )dxdt 

t=0 t=Q x=0 i=2 
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Or A J r = 


I ! r 

/ 

1=0 L 
'/ 1 


+ 


A 1 < (•'/:, <) 


ATr 


^(A7(;) 

Ox 


t = 0 x=Q 
'/ 1 

- / /AC 


1 w 

1=0 t ^) 

5AM 


,f - r n\ 1 1 

AT c K(x,t)A- 


dt, 


J x=0 


l (*<*•% 


dxdt 


(=0 i=0 
*/ 1 


A 


ar 

at 


+ 



f=0 x=0 
'/ 


a(CA) 

dt 


i 

dxdt — J [C XATcftLodx 

x=0 

h 

ATcdxdt + 2 / [(T- F) AT c ]x=i 

f=0 


+ 2 J [(T - F)AT c ] x=0 dt 
/.=() 

*/ l A/ — 1 

+ 2 J j Y'lT-Y] AT c 6{x - Xi)dxdt 


(2.85) 


i=0 x=0 


— n t=2 


From the sensitivity equation (2.51-2.52) we have 

at x = 0 and x = 1 




dx 

Hence the first, integrand of Equation (2.85) will vanish. Thus we have 
'/ i 


A J c 


= // A7V 
i=0 x=0 

+ I AT < 


o_ 

dx 


K(x,t) 


dX\ d{C A) 


ax 


+ 


at 


M — 1 

+ 2 £ P* - Y W - Xi) 

i = 2 


dxdt 


t-A) 

If 


OX 


K(x,t)g^ + 2(T-Y) 


dt 


Ji=0 


/ A7; - 


f .=0 
1 


A-(x, ( )|-2(7'-y) 
*/ 1 


dt 


J X=1 


J [CXATciUdx - j f AC 


x=0 


t = 0 x=0 


A 


or 

Ox 


dxdt 


( 2 . 86 ) 


Note that A Jc should not depend on A Tq. Therefore the integrands containing ATc are 
zero. Thus we have 
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ATM)— +2(T-Y)=[) 

for x = 0 

(2.88) 

K(x,t)^-‘2{T-Y)=0 

ax 

for :r = 1 

(2.89) 

A = 0 for t 

= t f 

(2.90) 

From the sensitivity equation (2.53) we see that ATc is zero at t = 0. 

thus reduced to 

Equation (2.86) is 

^ =-11 [ A ! 
t-t) x=() L 

ACdxdt 

(2.91) 

Comparing the Equations (2.79-2.82) and (2.87-2.90) we see that the adjoint equations 

for I< and C are similar. 

The perturbations A J K and A J c can be defined as 


t S 1 

AJk = J J J'x(x,t) AKdxdt 

l-{) x-() 

(2.92) 

'/ i 

A J c = f \ J' c {x,t)ACdxdt 

t= Ox-O 

(2.93) 

Comparing (2.83) and (2.92) we get the following expression for the gradient of the func- 
tional with respect to K: 

P t n 0X0T 

J k(x, t) = - — — 
ox ox 

(2.94) 

Comparing e(juat,ions (2.91) and (2.93) we get the gradient of the functional with respect 
to C: 

J (, (x, f ) A(x, t) 

(2.95) 


The gradient search procedure for obtaining the property functions K and C can be 
constructed as follows: 


1. Assume the properties K (:r, /.) and C{x, /.); In the absence of additional information, 
they can be treated as constants. 
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2. Generate the corresponding temperature field T(x,t), by solving Equations (2.42- 
2.45). 

3. Solve the adjoint problem, Equations (2.79-2.82), and obtain adjoint variable A (x,t). 

4. Compute J' k and J'c from Equations (2.94) and (2.95) 

5. Calculate the conjugate coefficients i/r and vc using Equations (2.59) and (2.60 ) 
respectively. 

6. Estimate the directions of descent P K and Pc from Equations (2.57) and (2.58) 
respectively. 

7. Solve the sensitivity problem, Equations (2.46- 2.49) and (2.50-2.53) to obtain sen- 
sitivity functions, 3ST K and A Tc- 

8. Compute step sizes (3k and (3c using Equations (2.73-2.74 ). 

9. Estimate K and C using Equations (2.55-2.56 ). 

10. Repeat the above calculation procedure until the discrepancy principle given by 
Equation (2.36) is satisfied. 

2.3.2 Flux- Temperature Boundary Conditions 

Direct Problem 


Consider tin? following transient inverse heat conduction problem. A slab of thickness L is 
initially at temperature f 0 . The boundary surface at x = 0 is subjected to constant heat 
flux qi while the other boundary at x = L is insulated. Then the GDE of the physical 
problem can be written as 



= f>C(T) 

dT(x, t) . - 

— y-L in 0 < X < L 

at 

(2.96) 

ox 

= Qi 

at x = 0 

(2.97) 

f 

= f t 

at x = L 

(2.98) 

T(x,t) 

= To 

for t — 0 

(2.99) 
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Assuming the following dimensionless quantifies 



X 

= V 

H' 

/: ■ A IT, 


z 

K r t 

Lqi 

r- C T- T ‘ 

c ~c , :■ 1 ~ % 


~ pC r l 2 ' ^ 

' KrV 


we obtain dimensionless form of e< 

{nations (2.96-2.99 ) as follows 


-( 

Ox \ 

v OX 

) = Cm° T ^ t] in0<x<l 

(2.100) 


-K { T) dT ^ t] - « 

at x = 0 

(2.101) 



T = Ti 

at x — \ 

(2.102) 



T = To 

for t = 0 

(2.103) 


Here the superscript. and subscript ‘r’ denote dimensional and referenced quantities, 
respectively. The sensitivity functions, adjoint functions and gradient equations for K 
and C are adopted as follows: 


Sensitivity Problem and Search Step Size 


The sensitivity function A T K is the change in estimated temperatures at different positions 
and times due to change A K in K(x,t). The sensitivity function AT# is given by: 

,d AT K (x,t) 




0 (a r-/ .X dT{x,t) 


= c- 


dx 


I\(x,l 


OATk (x, /.) 


A K(x,t) 


OT{x, t) 


Ox , " / Ox 

A1'k(x, t) = 0 for x = 1 


for x — 0 


A7’, v (:»:,/.) rr () for/, = () 

Sensitivity function A Tc is determined from the equation: 


for 0 <(2.104) 

(2.105) 

(2.106) 
(2.107) 


d ( K{x,t) ^^^ \ +AC(r,f)^M for 0 < £ < 1 (2.108) 

OX 


dx 


dx 


dx 


0ATc(x, t) 

Ox 


= 0 for x = 0 


(2.109) 


A T c {x, t) = 0 for x = 1 


( 2 . 110 ) 
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AT r {x, t) = 0 for t = 0 


( 2 . 111 ) 


The objective function can be written as before: 


r m 

J[I<(T),C(T)} = J[K{x,t),C(x,t))= / (2.112) 

tl o *= i 

where, 1 /,) is the temperature reading taken with a sensor placed at x = Xi at time t. 
M sensors are used to record the temperature history over time tf. Once again, all the 
sensors are assumed be uniformly placed. 

The following iterative process based on the conjugate gradient method is proposed 
to update K and C: 


K n+l (x,t) = K n {x,t) - P n K P%(x,t) 

C n+ \x,t) = C n {x,t)-^PS(x,t) 


(2.113) 

(2.114) 


The directions of descent for K and C are as follows: 


P2(x,t) = J' n c (x,t) + vZP2~\x,t) 


(2.115) 

(2. 11C) 


The sensitivity coefficients for K and C are given by 

/ f {J'lfdxdt 

, ,n f=0 i=0 

Vk ~ ~ ; 

f / (.I'^dxdt 

t -0 x=0 

f j ( J'lfdxdt 

n t ~ 0 x— 0 

c — tj , 

f / {.J'l^ydxdi 

t=0x=0 


(2.117) 


(2.118) 


Adjoint Problem and Gradient Equation 

To derive the adjoint problem for K(x, t), Equation (2.100) is multiplied by the Lagrange 
multiplier A(.r, t) and the resulting expression is integrated over the time and space do- 
mains. Then the result is added to the right-hand side of Equation (2.112). Hence the 
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esulting expression of the functional J[K(x,t),C(x,t)} becomes 

b i 


J[K(xJ),C(xJ.)} = I j X 


t-i) , T — 0 




dxdt 


r rn 

+ J O Ti{xi,t) -Yi( x i,i)f 


dt 


(2.119) 


t=o i=1 


The variation A J A is obtained by perturbing T by A T A in Equation (2.119), subtracting 
from the resulting expression the original Equation (2.119) and neglecting the higher-order 
terms. We thus obtain 


b i 


0 (... ., d(AT x ) \ a ( dT\ 3AT K 

~'dx at 


a W / a sr' ,) fc + £ 

/-()x=0 L \ / 


dxdt 


h i 


+2 j j 'jp L T-Y]AT K 6(x-x i )dxdt 


( 2 . 120 ) 


v 47 ; i 

t=0x=0 z -‘ 


Or A J K = 


<i r 
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t—0 

<1 
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/.=0 
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A K{x, t ) 


9(AT k ) 


9a; 
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0 i 


AAA' 


.9T 


dx 


dt 

J i=o t=0 x=0 
l b i 

dt 

1=0 i=0i=0 

b J 


J J K (x, t) 


OX 

dx 


9(A T k ) 


da: 


dxdt 


-//-II 


dxdt 


- f [CXAT K ]& Q dx+ J J AT K ^p-dxdt 

x=() f=0 1=0 

'/ '/ 

+ 2 / [(T-y)AT*] x=1 dt + 2 [ [(T-Y)AT K } x=0 dt 

t = o t=o 

b ' M—i 

+ 2 j J J2 i T ~ Y]AT k 5{x - x { )dxdt 


t=0 x=0 


Or A Jk- = 
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f=0 


A A' (a:, t) 


d(AT A ) 


da; 


d,t, -b 

J 1=0 (=0 x=0 


b i 

/ / A7V 




dxdt 
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' I r 
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t 0 


A7’„A'(:M)r 


OX 


Ox 


J x O 


/-:() 


f f r 'VTll 1 

aaa-A 

ax 


U i 


( i i 


,11- j j AK 

J.T -0 /”“() ,T“() 
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a\ar 

0.T C?X 


dxdt 


— J [C XATf{\ t L 0 dx + j j A7/< — — — dxdt + 2 J [(i — i ) Ai k]x=i 

x=0 t=0 x=0 t=0 

*/ *7 1 M-l 

+ 2 j' 1(7’ - iWKWorf* + 2 J j El 7 ’ - V']A7V<5(.r - Xijdzdt 
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i o 


f— 0 1=0 


Or AJ a - = 


'/ i 

II 


A Tn 


/,=() X=() 

+ J at k 

t= 0 
1 


<9:r 


+ ^ + 2 g [t - - *<> 


dxdt 


d\ 


K(x,t)- + 2(T-Y) 


*/ 

dt — J A T# 

J x=0 ^_o 


SA 


/<(x,t)— -2(T-F) 


dt 


X—l 


^ a 1 + a ^ 

ax ax 


J [C\AT K ] t t L 0 dx+ J \\K{x,t) ^=^ + AK 

=0 

l f 1 

/ / A* 


dt 


x=0 


aAar 

ax a:r 


dardi 


( 2 . 121 ) 


t=o i =0 

Prom the sensitivity problem at the two boundaries we get 


I<(x, 0 --^— + A/C^- = 0 (2.122) 

As stated earlier, A./^- is not a function A7 ’k. Therefore the integrands containing AT# 
should be zero. Hence 


a / , N aA\ a(CA) 
fc( A( *''b,) + « 

M-l 

+ 2 ^['T-y]<5(x-x i ) = 0 

i= 2 

(2.123) 

A ( x,/,)- + 2 ( T- 

T) = 0 at, x = 0 

(2.124) 

A = 0 

at x = 1 

(2.125) 

A = 0 

at t = tj 

(2.126) 
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’he Equation (2.121) reduces to 


A./;, 


// 


A K 


t = 0 i=0 


dXdT 
Ox Ox 


dxdt 


(2.127) 


ly definition 


A J, 


K 


>f 1 

J J J' K dx dt 


(2.128) 


(=0 i=0 


(2.129) 


Comparing the Equations (2.127) and (2.128) we get 

V - J^JL 
K dx dx 

similarly, considering the variation of C , we obtain the following expression for the gra- 
lient of the functional with respect to C: 

. dT 


J'c 


dx 


(2.130) 


The gradient search procedure for obtaining the property functions K and C can be 
constructed as follows: 


1. Assume the properties K(x,t) and C(.r, t)] In the absence of additional information, 
they can be treated as constants. 

2. Generate the corresponding temperature field T(x,t), by solving Equations (2.100- 
2.103). 

3. Solve the adjoint problem, Equations (2.123-2.126), and obtain adjoint variable 

X(x,t). 

4. Compute ,J' x and J'c from Equations (2.129) and (2.130) 

5. Calculate the conjugate coefficients u K and vq using Equations (2.117) and (2.118 
) respectively. 

6. Estimate the directions of descent I\ and Pc from Equations (2.115) and (2.1 1C) 
respectively. 

7. Solve tlie sensitivity problem, Equations (2.104- 2.107) and (2.108-2.111) to obtain 
sensitivity functions, A Tk and A Tq. 
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8. Compute step sizes ft]< and flc using Equations (2.73-2.74 ). 

9. Estimate' K and C using Equations (2.113-2.114 ). 

10. Repeat the above calculation procedure until the discrepancy principle given by 
Equation (2.36) is satisfied. 

In summary, the change from flux-flux to flux-temperature boundary conditions leads to 
the following difficulties: 

1. In case of flux-flux boundary conditions, the direct problem is always unsteady 
whereas it becomes steady a after certain time for flux-temperature boundary con- 
ditions. 

2. The adjoint function becomes zero at the boundary with the fixed temperature. 
This can affect the inverse solution. 


2.4 Coupled Equations without Source Term 


The inverse procedure for two inter-dependent space variables is developed in the present 
section. For definiteness, the source terms have been set to zero. Inversion of data with 
source terms is considered in Section 2.5. 


2.4.1 Steady State Problem 

Direct Problem 


Consider the following set of coupled equations for state variables Tj and T 2 properties 
K] and Ay. 


d_ 

dx 


k x {i\,r 2 ) 


'EM ] 

dx ) 


= 0 


(2.131) 


_d_ 

dx 


K 2 (T U T 2 ) 


m(x) \ 

dx ) 


= 0 


(2.132) 
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dl\ 

A i - - - q | 

ax 


at x = 0 


(2.133) 


= " 2 


at x = 0 


2’, (x) = T lr , T 2 (x) = T 2r at x = 1 


(2.134) 

(2.135) 


Inverse Problem 

For the inverse problem, A' 1 (T 1 ,T 2 ) and A 2 (Ti,T 2 ) are regarded as being unknown, but 
everything else in equations (2.131-2.132) is known. The experimental values of T\ and 
T 2 at some appropriate locations are considered available. The solution of the inverse 
problem is to be obtained in such a way that the following functional is minimized: 

M 

J[KuK 2 ) = J[I< l {x),K 2 {x)] = £[Ti(*i) - Fj^)] 2 

1=1 


M 

+ ^[T 2 (x { ) - Y 2 {xi)] 2 (2.136) 

t=i 

where Fi(;r,) and l 2 (x ? :) are experimental values of 1\ (xj) and 'i 2 (x<) respectively. 


Conjugate Gradient Method for Minimization 

The following iterative process based on the conjugate gradient method is used for esti- 
mation of A'i(x) and A' 2 (x) by minimizing the above functional J[K\ (x), A’ 2 (x)] 


i<r'a) = k?(x) - Kxiiia 

(2.137) 

KF'(x) = k; (*) - P'kiPU*) 

(2.138) 


where ft )\ , , are step sizes for A'i and A' 2 in going from nth to (n + l)th iteration. The 
directions of descent, P]\ x and PjJ, for A'! and A' 2 are given by 

ru*) = J" Kl (x) + <%,P£'(x) 


(2.139) 
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(2.140) 


These are the conjugation of the gradient directions J' n KX and at iteration n and the 
directions of descent ’ and P at iteration n — 1 for K\ and K 2 respectively. The 
conjugate coefficients are determined from 

fx = 0 X {J' n K\) 2 dx 


»k i = 


fx = 0 WliF** 


with v° K j 


(2.141) 


'l<2 


I 

£=0 

f ( J' n K2) 2 dx 

x=0 


with i/j <2 = 0 


The intermediate steps in the inverse procedure are summarized below: 


Sensitivity Problem for K t 


(2.142) 


d_ 

dx 

d_ 

dx 


K \{x) 

I<2(X) 


d(ATi) 


K 1 


dx 

d(AT 2 ) K1 


dx 

h 1 <ix 


-It 


(A T 2 ) 


K 1 


dx 


d_ 

dx 


A kM- 

dx 


= 0 


A / • dTl 

A/v, 

dx 


0 


= 0 


(A7j) ;a =0, (A T 2 )ki = 0 


at x — 1 


(2.143) 

(2.144) 

(2.145) 

(2.146) 

(2.147) 


Sensitivity Problem for K 2 


d 

dx 

d_ 

dx 


Ki(x) 

K 2 {x) 


d{AT\)x2 


dx 

d(AT 2 )K2 


dx 


dx 


(AT 2 ) a -2 


= 0 


+ 


d_ 

dx 


A K, 


dJ 2 
dx 


(AT\)i< 2 n t n 

—hi : = 0 at x = 0 


v r- d.T 2 

AA 2 — at a: = 0 

dx 


{ATi ) K 2 = 0, (AT 2 )k 2 = 0 at x = 1 


(2.148) 

(2.149) 

(2.150) 

(2.151) 

(2.152) 
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Step Sizes 


The functional ,J(K j , ) for iteration n 4- 1 is given by 

Af 

■7(/<? + ', a? +, ) = EPiW - .*■? - ffiA) - i'll? 


i=l 


Af 


+ EtW - n,r n K„i<Z - K,rh) - U] 


(2.153) 


?:■ ! 


where 7t " +1 and A'." 11 have been replaced by the expressions given by equations (2.137- 
2.138). Expanding (2.153) using Taylor expansion and neglecting higher-order terms we 
get 


M 

E 


M 


re, 


„ ar, 


+ E 

Or .7(A; i+1 , Ap 1 ) = 


T 2 (A?,A 2 n )-/^ 


l 5Ai 

9t 2 
9 Ai 


BT 

&K2 ( Pk2’qt7~ | — U 


: 9A, 


^2 ' ^K2 


n 

dl<2 


J i 
2 


-f 2 


A/ 


Eir.w./cjj-ffi.Am)!, -feA(T,r K2 -ii];+ 

?“1 

M 

E [r 2 (Xf, A-J, C”) - ffi, A(T 2 )" , - H" K2 A(T 2 r K2 - K,lJ dt (2.154) 

t=l 

where A JO = P)(-, and AA’ 2 = P}' <2 - Now the step sizes , and /?£ 2 in such a way that 
the functional given by the Equation (2.154) is minimized. Therefore 

<9,7” +1 


i 

9.7 n+1 

W 


= 0 


= 0 


(2.155) 

(2.156) 


K2 


From Equation (2.155) we get, 


M 

ffi, E(( Ar i)5n + (AT 2 ) 2 k ,]„ 
1=1 
M 


+ fix 2 E[(^r , i)yi(ATi)K2 + (^r 2 )Ki(AT 2 )^ 2 ] I=2 


t=i 


A/ 


= EK T i - FiKATi)*! + (T 2 - y 2 )(AT 2 ) jn ] I 


(2.157) 
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Equation (2.167) can be 

written as 



Gil &K\ + a l20K2 = h 

(2.158) 

Similarly, from (2.156) we have 



«21 0K\ + 0,22/3k2 — b 2 

(2.159) 

where 

M 

EK at i)ki + 

1=1 


Oil = 

(2.160) 

Oi2 = 

M 

+ (^32)ki(AT 2 )k2]x=21 

1=1 

M 

(2.161) 

bi = 

£[(T 1 - Yi)(ATi)ki + (T 2 - y 2 )(AT 2 ) K1 ]x=2i 

1=1 

(2.162) 

d2\ = 

M 

YX^T X ) KX {AT X ) K2 + ( AT 2 ) K 1 ( AT 2 ) K2]x=xi 

1=1 

(2.163) 

«22 = 

M 

EK at i)k 2 + (AT 2 )^ 2 ] I=Ii 

1=1 

(2.164) 

& 2 = 

M 

ym - >i)(a t x ) K2 + (t 2 - r 2 )(AT 2 ) K2 ]x=x« 

(2.165) 


i=i 


Adjoint Problem 


To derive the adjoint problem for K\(x) and K 2 (x) equations (2.131) and (2.132) are mul- 
tiplied respectively by the Lagrange multipliers A|(.t) and A 2 (.t) and the resulting expres- 
sion is integrated over the space domain. Then the result is added to the right-hand side 


of equation (2.136) to yield the following expression for the functional J[Ki(x), K 2 {x)\. 


J[K t(x) ,K M] = A, [;| (*,(.,^1) 

x=0 L \ / J x=0 L \ / 


x=0 

M 


dx 


M 


+ ~ ^it) 2 + T.(T 2l - Y 2 i) 2 (2.166) 

1=1 1=1 

The variation of the functional, A J k i due to change of parameter K\ to K\ -I- Aiifx is 
given by 


AJ 


Kl 




2=0 


d / d(AT K1 y 

Tx [ Kl{x) ^i^- 


d ( dTi 
+ di[ AK '-te 


dx 
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+ / A, ~ dx + 2 £ BT, - yi](AT,)/fi«5(a; - xjdx 

r L 0 L ’ \ *' / J x ” i^l 

/• 1 *" 

+ 2/ (2.1G7) 

Jx=0 7~i 


Or AJa 


A,A',(x) 




i i 


/ [ 

~n t- 


rfA x [d(ATi)Ki 


c£7 1 1 f \dX\ dT\l 

+ XxAKi— 1 ~ / A Ki ———— dx 

dx _ n J dx dx 

L X — 0 £~0 *” ’■* 

, f, , ^d(AT 2 )/ fl ] 1 f ,, , ,dA 2 [d(AT 2 ) KI 1 

+ XJU(X) ^— „ - / ~ 

L J i=0 x=0 L J 

+ 2[(T, - F 1 )(AT 1 )^ 1 ] I=1 + 2[(T 2 - F 2 )(AT 2 ) K1 ] I==1 
+ 2[(Ti - F 1 )(AT 1 ) ffl ]x=o + 2[(T 2 - F 2 )(AT 2 ) K1 ] I=0 

rl M-l 

+ 2/ ^[Tx-^KATj^ar-xOd® 

7x=0 j_2 
, A/-1 

+ 2 / E [T 2 - F 2 ](Ar 2 ) K1 6(x - X t )di 


Or AJ k1 


\ r,' t ^(ATOki ] 1 [ca-t\ t,- t \ 1 1 
AiAi(x) — — (ATi)j < iKi(x)-j— 

ax J x=o L ax J i=o 

r <i /„ , ,dA,M , . r, . „ 


/(A3-,)k, 


d ( T , , ,d\i\ , . a r *■ dTi 

_ A| ( X) _ *,+ A,AA,^ 


f ,\f/ fdA.dTil [ d(AT 2 )^i]' 

/ AAl rf7rf7 dx+ AgA2(l) dj 

W UulV li/J/ tlU./ A 

1=0 L J L J x=0 


(AT 2 ) /a A' 2 (x)- 


x=0 , 


/ {AT 2 )ki i[K^ 


+ 2[(T t - FXA'rO/nUi + 2[(r 2 - F 2 )(A7 2 )k,]x=i 
+ 2[(T] - \\){Al\) KX ] x=0 + 2[(7 2 - F 2 )(A7 1 2 ) a - 1 ] i= o 


rl A/-1 

+ 2 E^i- F 1 ](AT 1 )ki5(x - **)<& 

7x=o - =2 
rl Af-1 

+ 2 / Em- F 2 ](AT 2 )/a<5(x - x { )dx 

J x=0 • O 


(2.168) 


(2.169) 
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Now we have (A7 ’i)ki=(A7 2 )/(-i=0 at, x — 1. Thus Equation (2.169) becomes 
A .//vi = (A7’| ) Ki 


2(7’, -V',) 


+ 


l 

/(AT,) 

x~i) 

\ 

j (A7i) 


,T — 0 

M — 1 


+ I i Ai i 2 ^ K I 
r 0 


+ (AT 2 )/vi 


A' 


(lx 

+ I'm ->■>)«*-*<) 

rf(A 2 ) 

dx 


dx 


dx 


+ 2 (t 2 — y 2 ) 


+ 


+ 


A, (^(xj ^i d-A^f 


x=0 

1 


A 2 A' 2 (x) 


d(Ar 2 ) yi 

dx 


J x=0 T _, 


/ AA, 
x — 0 


J x=0 

dAj dT\ 
dx dx 


dx 


(2.170) 


We assume that the values of Lagrange multipliers A! and A 2 are zero at the two bound- 
aries. Since A J KX does not depend on the variations (AT) ) /n and (AT 2 );a, the integrands 
containing (A7) )/ <i and (A7 2 )/<i are zero. We thus have 


dx 

d_ 

dx 


M — l 


^A'i(x)^-j + 5 Z ( r i ” - *i) = 0 


,/ / \ M — l 

Ki~- + 2(Ti — Fi) = 0, at x = 0 
ox 

I<- 2 ~ + 2(7 2 - V 2 ) = 0, at x = 0 
dx 


(2.171) 

(2.172) 

(2.173) 

(2.174) 

(2.175) 


Ai = A 2 = 0, at x = 1 
Equations (2.171-2.175) can be solved directly to give the adjoint functions Ai and A 2 


Gradient Equations 


Using Equations (2.171-2.175) Equation (2.170) can be reduced to 

i 


AJj<- 


K 1 


- / AA > 
x=0 


dXi d'l\ 
dx dx 


dx 


(2.176) 
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From the definition of AJk\ we have 

1 

A J KX = J J 1 K)(:r.) A K\dx (2.177) 

x~() 

and comparing equations (2.176) and (2.177) we get 

J'k i = i (2.178) 

ax ax 

Similarly, considering the variation of the objective function due to variation of K 2 (x) one 
can get the following gradient equation 

f K2 = -A d Jl (2.179) 

ax ax 

The inverse technique for reconstructing K\ (7j , T 2 ) and K\ (Ti, T 2 ) can be cast in the 
following algorithmic form: 


1. Assume the properties K\{x) and A' 2 (x); In the absence of additional information, 
they can be treated as constants. 

2. Generate the corresponding state variables r l\(x) and T 2 (x) by solving Equations 
(2.131-2.135). 

3. Solve the adjoint problem, Equations (2.171-2.175), and obtain adjoint variables 
Ai (x) and A 2 (:r) . 

4. Compute J'/<i and J' k 2 from Equations (2.178) and (2.179). 

5. Calculate the conjugate coefficients vki and vk 2 using Equations (2.141) and (2.142 
) respectively. 

6. Estimate the directions of descent Pki and Pk 2 from Equations (2.139) and (2.140) 
respectively. 

7. Solve the sensitivity problem, Equations (2.143- 2.152) to obtain sensitivity func- 
tions, A T 2 k\i ATi, < 2 and A T 2 ;< 2 - 

8. Compute step sizes /Aa and (3k2 using Equations (2.158-2.159 ). 

9. Estimate K\ and K 2 using Equations (2.137-2.138 ). 

10. Repeat the above calculation procedure until the discrepancy principle given by 
Equation (2.36) is satisfied. 
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2.4.2 Transient Problem 

Direct Problem 


Consider the system of equations 




dT x (x, t) 

-c ( r,,r 2 ) m 

(2.180) 

|( A , (Tl , T2) EMM) 

nr r\ dT ^ x,t ) 

i =C(T i ’ T2) dt 

(2.181) 


= qi at x — 0 

(2.182) 

„K 2{TuT2) d J^ 

= q 2 at x = 0 

(2.183) 

II 

at x = 1 

(2.184) 

T 2 (x, t) = T 2 l 

at x — 1 

(2.185) 

II 

H 

for t — 0 

(2.186) 

T\{x, t) = l\o 

for t = 0 

(2.187) 


with Ki, K 2 and C as the undetermined parameters. 

Inverse Problem 

For the inverse problem, K\ (T\,T 2 ), K 2 (Ti,T 2 ) and C(l\, T 2 ) are regarded as being un- 
known, other quantities in Equations (2.180-2.187) being known. The experimental values 
of T] and 7 2 at some appropriate locations are considered available. The solution of the 
inverse problem is to be obtained in such a way that the following functional is minimized: 

l f M 

J[K u I< 2 ,C] = J[K\(x, t), K 2 (x,t), C(x,t)] = j - Y x {x h t)\ 2 dt 


t = o i=1 


7 M 


+ [ 

do i=1 


dt 


(2.188) 


where Y x (xi,t) and Y 2 (xi,t) are experimental values of T x (xi, t) and T 2 (x t , t) respectively. 
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Conjugate Gradient Method for Minimization 

The following iterative process based on the conjugate gradient method is used for estima- 
tion of A'i (:r, /,), K 2 (x, /,) and C(x, t.) by minimizing the above functional J[Ki(x, t),K 2 (x, t),C(x , £)] 



(2.189) 

K +l (x,t) = A 2 "(x,i) - 

(2.190) 

C" +I (.x,t) = 0 , t) - F c PZ(x,t) 

(2.191) 

where t) , fi)\ 2 {x, t) an d (3£(x,t) are step sizes for K\, K 2 and C in going from the 

nth to (n + 1 )th iteration. The directions of descent, P]q, Pj^, Pq , for K 1 , K 2 and C are 

given by 


CjritM) = 

(2.192) 

PJ 2 (x,() = ./™ 2 (x,() + ^ 2 PE'(x,i) 

(2.193) 

PS(x,t) = J'Hx,t) + ^Pg-'(x,t) 

(2.194) 

These are the conjugation of the gradient directions J'ki an d J'c iteration n 

and the directions of descent P£ 7 1 , P) l <2 1 and P}P X at iteration n — 1 for K\, K 2 and C 

respectively. The conjugate coeflicients are determined from 


v , ; sum with „o 

Al sUJ'h'rdx M 

(2.195) 

- J*=o( J "K2) 2 ' Ix _ j ( | 0 _ 

Ka ~ " lth ‘ V2 “° 

(2.19G) 

,, _ Willi wf. = 0 

/.UO'c ') 2 <& 

(2.197) 


Sensitivity Problem and Step Size 

There are three unknown parameters K i, K 2 and C which are to be determined. In order 
to derive the sensitivity problem for each unknown, we should perturb the unknowns are 
perturbed one at a time. 
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Sensitivity Problem for K\ 


It, is assumed that when K\(x,t) undergoes a variation AK\, l\(x,t) and 7 2 (x,f) are 
perturbed b.v 7’, + A 7’, and 7 2 + AI 2 respectively. Then replacing K\ by K + K \ , T x 
by T i + (A7 'i)ki and 7’ 2 by 7 2 + (A7 2 )/<, in the direct problem, and subtracting from 
the resulting expressions and neglecting the higher-order terms, the following sensitivity 
functions (A7’ 1 )/ V | and (A7 2 )/vi are obtained. 
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Ox 


I< 


d(A r,) K1 


Ox 
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tlx 


A/i i 


dTi 
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d(A T,)*, 


Ox 


3(AT 2 ) Kl 
; ' 2 ax ' 
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Ox 
d(AT 2 ) Kl 


Ot 


K\ + AKl °li = 0 

ox ox 


Ot 

at x = 0 




(AT,),,, =1) 
{AT\)k\ = 0 


Ox 

and 

and 


at x — 0 
(A7 2 )/n =0 at x = 1 

(AT 2 )ki — 0 at t = 0 


(2.198) 

(2.199) 

( 2 . 200 ) 

( 2 . 201 ) 

( 2 . 202 ) 

(2.203) 


Sensitivity Problem for /i 2 


Following an identical procedure in deriving the sensitivity functions for K\ we get the 
sensitivity problem for K 2 as follows 


0_ 
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k 0{AT x ) K2 


dx 


= C 


0(A1\) K2 


Ot 
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y(A7 2 )/< 2 

iY‘> ^ 
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+ 77 — 

r A ,, d7’ 2 
AA 2 ~ 

ax 

L j 

Ox 

9 x_ 
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,<7(A7 2 ) A2 


Ot 


3(A7-,) K2 =0 

Ox 


at x = 0 


a(Aii)« Mi = 0 

ox ox 


at x = 0 

(AT,)/, - 2 = 0 and (AT 2 ) a - 2 = 0 at x = l 

(AT 1 ) /f2 = 0 and (A7 2 ) K2 = 0 at t = 0 


(2.204) 

(2.205) 

(2.206) 

(2.207) 

(2.208) 
(2.209) 
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Sensitivity Problem for C 

The sensitivity functions (A7’i)r.' and (AT 2 )c derived as 

0 \ v 5(A'A) c 1 ^d(A T x )c , 

01 hl ~dT~ = c ~~¥~ + l 


(Ar,)c: 


r, a(AT,)c 1 

' & j 

,=c a ^ + AC%' 
dt dt 

(2.210) 

a(Ar 2 )d 

a* 

d(AT 2 ) c , r dJ\ 

C dt ° dt 

(2.211) 

a(AT,)f 
A ' ax 

= 0 at x = 0 

(2.212) 

K B(AT 2 ) c 
h 2 Ox 

= 0 at x = 0 

(2.213) 

0 and 

(A T 2 ) c = 0 at x = 1 

(2.214) 

= 0 and 

(AT 2 ) c = 0 at t = 0 

(2.215) 


Step Sizes 

The functional J(A'" +1 , A'£ +1 ,C n+1 ) for iteration n+ 1 is given by 

l l u 

J(jrr + 1 ,A-J +, ,C" +1 ) = J Y;iTi(K? - 0 n K ,P^,K? - K 2 PZ 2 ,C" - ftPS) -Ytfdt 


/ M 

+ / ZlT0Kl‘-0hPZi,I<Z-PhPh,C"-0£PS)-Yi)ldt (2.216) 

do i=l 

where A'" 41 , K 2 l+t and C n+1 have been replaced by the expressions given by equations 
(2.189-2.191). Expanding (2.216) using Taylor expansion and neglecting higher-order 
terms we get, 

J(A'" +1 , A'.’,' 1 ',6'"") = 

/ E h(AT. JtJ.c*) - fit, (p^§-) - ffi 2 (pr 2 §-) - 01 (c5§) - Vif* 


/ E r.tAT.A-y.cn-ffl, 

do <=> L 

/EfctAr.Aj.c")-®, (pt,; 


ph Phgj - ) - ® l ps- 


Y 2 dt 
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Or ic;' \c n ") - 


r A1 

J E I'/’. (AT, AT,<7") - /^,A(7’,)k, - /^ 2 A(7’,)" 2 - ^A(T,)?; - H]? rft + 


<T, AT, C") - /?£, A (T 2 )” , - fie 2 A(T 2 ) n K2 - ft n c A {T 2 ) n c - Vi]f ^2.217) 


where AA'i = 7^,, A7 y 2 = 7^ 2 and AC = 7^. Now the step sizes, /7^-j, /j£ 2 and Pc, 
such a way that the functional given by the Equation (2.217) is minimized. Therefore 


dJ n+1 


(2.218) 


dJ n+l 


(2.219) 


dJ n+x 


( 2 . 220 ) 


From Equation (2.218) we get 


f M 

ft, ] E [(AT.^. + IAJPI.U^I 


f M 

+ pK2 J J2[(AT l )Kl{m)K2 + (A72) K1 (A7: 


2J/<2 J*=xit 


tin + (AT 2 ) c (Ar 2 ) K\]x=xidt 


An t=l 


■ M 

J EK 7 ’i - 1 ^(A'i’OK, + ('7 2 - l 2 )(A'/: 2 ) K! ] I=lif 7/. 


( 2 . 221 ) 


Equation (2.221) can be written as 


( h\ftio + d\2^K2 + a va ft a — C\ 


( 2 . 222 ) 


Similarly, from (2.219) and (2.220) we have 


&21 @Kl + 022 ft K2 + 0 2 $/7c — C 2 


(2.223) 
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0.31 lh< 1 + 0 , 32 /^ 2 + 0 , 33 /ic — C 3 

(2.224) 

where 



O-ll 

'/ M 

= / B(^i) 2 ki + (AT 2 ) 2 K1 ] x=li dt 

£=0 * =1 

(2.225) 

n i2 

*/ M 

= [ f^[(AT,) KX (AT x ), <2 + (AT 2 ) Kl (AT 2 ) K2 ] x=xi dt 

t=0 

(2.226) 

Oi3 

*/ M 

= 1 'Z[(AT i ) K1 (AT 1 ) c + (AT 2 ) K1 (AT 2 )c] x =xidt 

d 0 i=1 

(2.227) 

Cl 

7 M 

= / Elm - >i)(ATi)ki + (T 2 - F 2 )(AT 2 ) K i]x=xidt 

t=o ’= l 

(2.228) 

021 

7 M 

= J EKATijiaCAT:)^ + (Ar^^Ar,)*,]*^ 

t=o i=1 

(2.229) 

0-22 

7 Af 

= / EK A:r >)K2 + (&T 2 )U*=,idt 

1=0 i=L 

(2.230) 

«23 

= / EKATO^IATOc + (AliMATiM^dt 
f=0 i=l 

(2.231) 

02 

7 A/ 

= / Elm - r,)(AT ,) K2 + (T 2 - y 2 )(AT 2 ) K 2 W^ 

/' () 1 

(2.232) 

0.3 1 

7 A/ 

= / / ,[( AT) ) k 1 ( A? 1 )g + ( AT 2 )k\ (A -l 2 )c]:r=3;x<it 

t =0 i=1 

(2.233) 

0.32 

t S M 

= [ El(AT I ) K 2 (AT 1 ) c + (AT 2 ) K 2 (Ar 2 ) r; ] x=li dt 

/ 0 <=> 

(2.234) 

033 

*7 

= [ U(^)b + (AT 2 )l] x=xi dt 

t L 0 i=l 

(2.235) 

0.3 

7 A/ 

= / Elm - r,)(AT,) c + (T 2 - Y 2 )(AT 2 )c} x = xi dt 

do i=» 

(2.236) 


Solving Equations (2.222-2.224) simultaneously we get step sizes Pk\, Pk 2 and fi c . 
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Adjoint Problem and Gradient Equations 

To derive the adjoint problem, Equation (2.180) is multiplied by the adjoint function 
\i(x, t) and Equation (2.181) is multiplied by X 2 (x,t). The resulting expression is inte- 
grated over the space and time domains. Then the result is added to the right hand side 
of Equation (2.188). Then the objective function becomes 


t= o i=1 


J[K U I< 2 ,C] = J[I< 1 (x,t),K 2 (x,t),C(x,t)} = I fyTiM-Yiixirffdt 
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+ / / A] j — | h ] — | — C 

Mil'll 1 


dx 
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l ! 1 


+ 




0 x 


j j 

1=0 x=0 


dx 


dt 


dxdt 


dxdt 


(2.237) 


The variation A Jk\ is obtained by perturbing Tj by (ATi)ki and T 2 by (AT 2 )k 2 in 
Equation (2.237). Then the original Equation (2.237) is subtracted from the resulting 
equation. After neglecting the higher-order terms we get the following expression of the 
variation A J K \\ 

b a/ b 
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(2.238) 
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I '7 » 
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(2.239) 


Or A J Kl = 
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t— 0 x=0 
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(2.240) 
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Now the integrands cont aining ( Al \ ) /<• j and (AT])ki should he zero. We thus have 


0 (,. OXA , 0(C A,) A/ -' 


ihr V ' 1 ih. 


4* 


Dt 


+ 2 £ [T, - - .t { ) = 0 


i-2 


d\ 

K\ -j-*- + 2 (Tj - li) = 0 for x = 0 

ox 


Aj = 0 for x = 1 


(2.241) 

(2.242) 

(2.243) 


Ai = 0 for t = t< 




BX 

I< 2 ~ + 2(72 - V 2 ) = 0 for z = 0 
ox 


X 2 = 0 for x = 1 


(2.244) 

(2.245) 

(2.246) 

(2.247) 


A 2 = 0 for t = tf 


b i 

A Jki = ~ j j A K x 


From the definition of A Jki 


t — 0 .T — 0 


If 1 


dXi dTi 

dx dx 


dxdt 


(2.248) 


(2.249) 


A J/<i = f f J' k\ A A i dxdt 

t= 0i=0 

Comparing Equations (2.249-2.250) the following gradient equation is obtained 

„ , dX,dT x 

dx dx 

Considering the variation of K 2 one can show the following relationship 

ox 2 or 2 


(2.250) 


J K2(-T)^) 


dx d x 


(2.251) 


(2.252) 



2.4 Coupled Equations without Source Term 


54 


For deriving the expression for the gradient of the functional with respect to C, we need 
to apply perturbation principle in Equation (2.237). It is assumed that when C(x,t) 
undergoes a variation A Tj{x,t.) and /.) are perturbed by T\ + (ATi)c and 

7 2 + {AT 2 )c- Then replacing C by C -f AC, T\ by 'i\ + (A'i’i)c and T 2 by T 2 4- (AT 2 )c 
in Equation (2.237), subtracting Equation (2.237) and neglecting the higher-order terms, 
we get the variation of the functional as follows: 

b a/ b m 

AJ C = 2 j '£{l\(x i ,t)-Y l (x i ,t)}(AT l ) c dt + 2 J jT[T 2 fc, t) - Y 2 (x u t)](AT 2 ) c dt 


1-0 


i~ 1 


t—Q 


i=\ 


* 1 / 


+ 


/ - () x=Q 
'/ 1 

// 

t=0 i=0 


A, 


a („d{AT,)c\ 


dx 


K . 


9a: / 


_ Ar??k _ 

c at c at 


dxdt 


dxdt 


(2.253) 


Or A = 


+ 


2 J [Ti(xi,t) — Yi(xi,t)](ATi)cdt + 2 J [T 2 (x u t) - Y 2 (x u t)}(AT 2 ) c dt 

t = 0 t=0 

b b 

2 J [T x {x M ,t) - F x (x m , t)](ATi) c di + 2 j [T 2 (x M ,t) - Y 2 (x M , t)]{AT 2 ) c dt 


t - 0 

b l a / i 


+ 2 / f Y, [Ti (z, <) - 1 1(®, t)](A7\) c J(:E - *<)<&** 

(=0i=0 i=2 


'/ r */-« 

+ 2 j J £[T 2 (x,t)-V 2 (a;,i)](A T 2 )c6(x - Xi)dxdt 


f.=0 x—Q 7—2 

'/ i 


* / / Ml 


+ 


/. ();/: () 

'/ 1 

//*• 

/.- () x =0 

t f i 


o_ ( K 3(Ar 2 )e A _ ^, 0(Ar 2 )c 


9a; 


c- 


2 9a: j 

/ h c { x ^ +x ^) didt 

=0 x=0 x 7 


dt 


dxdt 


dxdt 


(2.254) 


Or AJc = 


'/ t 


/ / (A'A)c 


£=0 £=0 


y bC v 


9a: 


9a: 


t=2 


at 


dxdt 
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(A7 2 ) r 


0 (,, 0X 3 \ , 3(CA 2 ) 

01 [ h2 ~te r~~df~ 


+ (?2 ~ ^ 2 ) 5 (x - Xi) dx dt 


(Ar, M / ' , 7^ _2(ri_r,) ) dl 


(AT 2 ) c [K 2 -^-2(T 2 -Y 2 )\\ dt 


d(A T x )c 


dt— I A? A’ 


d(AT 2 ) c 


l 1 

- I \\iC(Al\) c }t f dx - I [\ 2 C(AT 2 ) c }tfdx 


AC ( Ai f + 


(2.255) 


Since AJ C floes not depend on (AT x )c and (AT 2 )c, the integrands containing (ATi)c or 
(AT 2 )c become zero. We thus liave 


0 (w. n h) + ?(£M 

<7x l 2 dx j dt 


+ x- Xi )=Q 

i = 2 

Af — T 

+ 5] ('ia - i: 2 )J(x - x<) = 0 


A i — 2{T\ — 1 1 ) = 0 at x = 0 

ax 

d\ 

I< 2 ~ - 2(T 2 - 1 2 ) = 0 at x = 0 

ax 

Aj = A 2 = 0 at x = 1 
Ai — A 2 = 0 at t = tf 


(2.256) 


Equation (2.255) now simplifies to 


t = 0 X—Q 




Prom the definition A Jr 


(2.257) 


AJc = J j J'c&Cdxdt 


(2.258) 
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Comparing Equations (2.257) and (2.258) the following gradient equation is obtained: 

(2.259) 


+ a,-' 


<)t. 


on 

<H 


The inverse technique for reconstructing Ki(Tx,T 2 ), K 2 (l\ , n) and C(2\, T 2 ) can be cast 
in the following algorithmic form: 


1. Assume the properties K\(x,t), Ii 2 (x,t) and C(x,t); In the absence of additional 
information, they can be treated as constants. 

2. Generate the corresponding state variables T\(x, t) and T 2 (x, t ) by solving Equations 
(2.180-2.187). 

3. Solve the adjoint problem, Equations (2.241-2.248), and obtain adjoint variables 
A i (:/: ) and A 2 (:r) . 

4. Compute J'ki, J' K 2 and J'c from Equations (2.251), (2.252) and (2.259) respec- 
tively. 

5. Calculate the conjugate coefficients iqn, V K2 and vc form Equations (2.195), (2.196) 
and (2.197) respectively. 

6. Estimate the directions of descent Pki, Pk 2 and Pc from Equations (2.192), (2.193) 
and (2.194) respectively. 

7. Solve the sensitivity problem, Equations (2.198- 2.215) to obtain sensitivity func- 
tions, (A7)) K1 , (AT 2 ) ku (A2\) k 2 , (AT 2 )k 2 , (A l\) c and (A7 \) c . 

8. Compute step sizes Ok i, 0k2 and flc using Equations (2.222-2.224 ). 

9. Estimate A',, A' 2 and C using Equations (2.189-2.191 ). 

10. Repeat the above calculation procedure until the discrepancy principle given by 
Equation (2.36) is satisfied. 


2.5 Coupled Equations with Source Terms 


The formulation of Section 2.4 is now extended to include source terms. The specilic 
choice of the source terms is motivated by the study of multi-phase flow in porous media. 
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2.5.4 Sensitivity Analysis 

Sensitivity Problem for A'| 


0 

dx 


I< 




dx 


d_ 

dx 


A K x 


dT x 


dx 


= c cH£T } )k 1 _ c d(AT 2 ) K1 


dt 


dt 


d_ 

Ox 


I< 


d(AT 2 ) K 


c 1 


dx 


c d(A T 2 )ki c d(A Ti) K i 


dt 


dt 


+ AK l ~ = 0 at x = 0 
ax ax 


(2.268) 

(2.269) 

(2.270) 


#2 


L J I\ I 

dx 

= 0 at 

x = 0 


and 

(A T 2 ) K1 

= 0 at 

x = 1 

and 

(AT 2 )k\ 

= 0 at 

t = 0 


(2.271) 

(2.272) 

(2.273) 


Sensitivity Problem for A ' 2 


a_ 

<9x 


^ d(AT l ) K2 


dx 


_ ^, d(A Tj)/<- 2 _ ( j9{AT 2 )k2 


dt 


dt, 


(2.274) 


c)x 


d(AT 2 ) K2 

iV o ^ 

d 

+ ~ — 

[a r- OTal 

AR 2 —A 

dx 

dx 

OX 


= c d{AT 2 ) K2 c d{AT,) K2 


dt 


„ 0 ( M \) K2 n 

A i = 0 at x = 0 

t/.r 


+ at * = 0 


dt 


dx 


dx 


(2.275) 

(2.276) 

(2.277) 


(AT\)k2 — 0 

and 

(AT 2 )k2 — 0 

at x = 1 

(2.278) 

(ATi) K2 = 0 

and 

{AT 2 )k2 = 0 

at t = 0 

(2.279) 
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Sensitivity Problem for C 


0 

dx 


h ' ~o-r~ 


c ^k + AC— - C°-^i _ AC m 


dt 


dt 


o_ 

dx 


0(A T 2 ) c 

h2 ~d^~ 


c a(A7’ 2 ) c , ^ g(A7^ 


dt 


+ C- 


c 


C 


Of. 


0(M\) C 


c?(A T\)c 

Ki v ~ = 0 


K, 


dx 

d(AT 2 ) c 

dx 


= 0 


dt dt 

at x = 0 

at x = 0 


- AC 


dt. 

dTi 


(A7j)c = 0 

and 

(A T 2 ) c = 0 

at x = 1 

(ATi)c = 0 

and 

(A T 2 ) c = 0 

at t = 0 


dt 


(2.280) 

(2.281) 

(2.282) 

(2.283) 

(2.284) 

(2.285) 


2.5.5 Step Sizes 

Equations (2.222-2.224) are applied to calculate step sizes, (3 K i, Pk 2 and (3c ■ 


2.5.6 Adjoint Problem and Gradient Equations 


To derive the adjoint problem, Equation (2.260) is multiplied by the adjoint function, Ai 
and Equation (2.261) is multiplied by the adjoint function A 2 and the resulting expression 
is integrated over the corresponding space and time domains. Then the result is added to 
the functional (2.188). We thus have the following expression of the functional: 


J[K u K 2 ,C] = .J[K l {x,t),K 2 {x,t),C{x,t)} = 


7 M 


t = 0 i=1 


J Y,[Ti{xi,t) -Y\(xi,t)] 2 dt+ j J2{T 2 (xi,t) - Y 2 (xi,t)] 

=0 *'~ 1 

f f 1 

//a, 


7 1 

+ 

t= 0 1=0 

7 1 

+ J J h 

t=0 x=0 


d ( v dTi(x,t)\ ry dTi(x, t) t n dT 2 (x,t) 
dx V 1 dx ° dt +c . dt 


'dt 


dxdt 


0 or 2 (x,t)\ , „ar,tM) 

ai [ h, ~ar~ r c ~d^ + bt~ 


dxdt . (2.286) 
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The variation AJki is obtained by perturbing 1\ by (ATi)ki and T 2 by ( AT 2 )k 2 in 
Equation (2.286). 'I’lien the original Ecpiation (2.286) is subtracted from the resulting 
equation. After neglecting the higher-order terms we gel, the following expression of the 
variation A./ki: 


AJki - 

1 A / 


+ 


2 j Ep\( 3 : i ,f.)-r 1 (x < ,t))(Ar 1 W< + 2 J 't\T 2 (x i ,t) -Y 2 (x u t)}(AT 2 ) Kl dt 

t - 0 j '~ ] t ~ 0 2=1 

l f 1 

//*■ 




AA'i 


f=0 x— 0 

l f 1 

+ J j A, 

f-0 x“() 

Or AJki — 


aiy 


a(AT 1 )Ki +c a(AT 2 )Ki 


dxdt 


d f 9(AT 2 )ki 

n 2 - 


9.x 


<9x 


_ c a(Ar 2 )Ki + ^d(ATi) K i 


/ 

t=o 

*/ 

/ 


1 \ SCAT,)^ 11 

Aj ill r 

OX 

„ C/Ai 


at 


f / i 


dt 


dxdt 


(2.287) 


dt + j J (^7’i)ki 


^Ok,*' 

- I ) AK 


f "0 

h 1 


1=0 t=0 x=0 

1 t l 


A f A- 
ax v 1 ax 


dxdt 


f- 0 x=0 
i 


dx 

dXi dT x 


dx dx 


M+J 
• T=0 f-o 

dxdt 

f f i 


fit 


Ji=0 


J \CMAT l ) Kl ] t t L 0 dx+ f j (ATOki ^^-dxdt 

x=() f— 0 x=0 


t— 0 x=0 

l f 1 


+ j [CX x {AT 2 ) m \\L 0 dx- j j (AT 2 )k, 


x— () 

1 r 


/ 

-0 

b 

- / 


f-0 
'/ r 


f=() 
l I 1 


X ^ 9(AT 2 )/C1 

^2^2 ~ 

ox 

(AT 2 ), fl K 2 A 


t=0 x=0 

1 b 1 

(it + 

Jx=0 t=0x=0 

1 1 


3(CAx) 

at 


dxdt 


— (l< ^ 

a.x l 2 dx 


I J (AT,). 

=0 x=0 
1 

dt- J [CX 2 (AT 2 ) Kl f t L 0 dx 


dxdt 


+ 


2)K\' 


d ^Mdxdt+ j [CX 2 (AT i ) K{ } t t L 0 dx 

X— 0 


£=0 x=0 

1 


f-0 X~() 
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+ 


+ 


+ 


'/ i. ‘)(CX 1 >! 

/ / + 2 y [(7’j — 1 |)(A7’i)Ki]x=ro 

I - II x 0 t = 0 

2 j K'i’, - Yi){ATi) Kl ] g=l + 2 

t=0 


/ M 

j J £['^i ~ ~ Xi){A'l\) Kx dxdt 


t= 0i=0 i=2 


t f t f 

2 f[(T 2 -Y 2 )(AT 2 ) KX ] x=0 + 2 j [(T 2 - Y 2 )(AT 2 ) KX } x=l 


t = 0 


f=0 


A-f 

2 / / - >2]r5(;r - x^AT^K.dxdt 

(=0 i=0 i=2 


(2.288) 


Or A,7/< = 
h i 


+ 


/ / ( Ar '> 

=0 x=0 

l f j 

/ (AT,) kx 


A'l 


A + - a J^l + 2 E ET. - V.W. - «,) 

i=2 


9a: V 9a: 


dt 


dt 


dxdt 


t= o 
i 




dt — f (ATi)ki 

x=0 t=o 




dt 


J X=1 


J [C\ ] (AT l ) l<l } l / = odx+ j \ x 
=0 

/ / (AT.) 


.T-~0 
l f \ 

+ / / 

(=0 i=0 

+ 

(=0 
tf 1 


a:r a:r 


dt 


x-Q 


°(k&) + - °-m + 2 eb - w. - «,)■ 


/ (Ar a ) K , 

' = 0 

l S 1 

/ / A *' 


Ox V * dx 


I<2^ + 2(r 2 - r 2 ) 


dt 


dt 


i = 2 


dxdt 


tf 


dt — J (AT 2 )k\ 

j x=0 t _o 


^_ 2( r 2 _y 2) 


dt 


x=l 


d\ x OTi 


0 „ . dxdt ( 2 . 289 ) 

I ax ax 

4=0 x=0 L J 

Now, the integrands containing (AT) ) ;<i and (AT 2 )ki should he zero. Thus we have the 
following adjoint equations: 

1 (*' A) + ^ ^ + 2 |V. - -*-)“» ( 2 . 290 ) 


a ^ + accM _ ^(cao + 2 ^ _ y . „ (i _ I() = 0 


9 a; V 9 a: 


9i 


dt 


( 2 . 291 ) 


i=2 


r)A 

A'j V- + 2 (Ti - Yi) = 0 at x = 0 
ao: 


( 2 . 292 ) 
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dX 2 
1 Ox 


A I — A-i — 0 


0 at x — 0 

(2.293) 

at :r = 1 

(2.294) 

for t = tf 

(2.295) 


OXi &T X 


dx dx 


dxdt 


(2.296) 


Equation (2.289) becomes 

b i 

AJ/n = — j J AAi 

l=0i=0 

From the definition of A J^\ we have 

b i 

A J Kl = j J J'ki^.OA/v, dx dt (2.297) 

l=Oi=(l 

Comparing Equations (2.296) and (2.297) the following adjoint equation is obtained. 

3A, dT x 


■J'kx = 


dx dx 


(2.298) 


The gradient of the functional with respect to Ji 2 is 

d\ 2 dT 2 . 

(2 ‘ 299) 

Let us derive the expression of J'c ■ For deriving the expression for gradient of the 
functional with respect to C, we need to apply perturbation principle in Equation (2.286). 
It is assumed that when C(x,t) undergoes a variation AC{x,t), T\(x,t) and T 2 (x,t) 
are perturbed by T\ + (ATi)c and T 2 + ( AT 2 )c ■ Then replacing C by C + AC, 7\ by 
T\ + (AT) )(7 and T 2 by T 2 -\-(AT 2 )c in Equation (2.286), and subtracting from the resulting 
expressions Equation (2.286) and neglecting the higher-order terms, we get the variation 
of the functional as follows: 

b M b hi 

Ada = 2 J^[T l {x u t.)-Y l {TiM*T l )cdt + 2 J ~ Y*M](AT 2 ) e dt 


t ,= 0 


i~ 1 


=n 


t = 0 


t f i 


Ih 

t - 0 T-0 

'/ i 

// 


+ 

/-() .T-0 

'/ I 

+ II X2 

t = 0 x=Q 


a („a{AT,)c\ r d(A T,)c a „9T 2 8(AT 2 )c 

oi \ K '~ a^) - AC ~af " C ~W~ + AC ~m + c at 


j>_ (,, ajA^k ' 

dx l 2 dx 


r a(AT 2 ) <; m\ r 0(M\)c 
AC - aT - c — m - + AC -aT + c — ar - 


dxdt 


dxdt 
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Or A J c = 


if if 

2 I |7i(.n,/) V| (.r t , 1)\(A'1\ )cdt + 2 f |7’ 2 (x,,t) - V 2 (x, , t)]( A7 2 Mt 

( o t=o 


+ 2 J [7\{x m J,) - Y\(xm, 

t = 0 f =0 

'/ !- A'-l 

+ 2 J J [Ti( x ’ t) - i ' r 1 (^, t )](AT 1 )c( 5 (x - Xi)dxdt 

t = 0 x=0 2=2 

*7 1 A?-l 

+ 2 / / [72(jc, «) - r 2 0M)](AT 2 )c<5(x - x { )dxdt 


2 {xM,t) - Y 2 {x M ,t)](AT 2 ) c dt 


i=0i=0 
'/ 1 


+ 


+ 


+ 


//A, 

rr() X~() 

'/ i 

If 


0 ( R d(ATi) c \ „0(A T,) c , „0(AT 2 )cl 


[9.x \ 1 9.x y 


it f K2 8(AT 2)c 


C 


9£ 


+ C- 


9£ 


x=() 

'/ l 

// 

t=0 x=0 


AC 


9.x 
(A, - A 2 ) 


dx 

djTjj-T]) 

dt 


c a(AT 2 ) c | c a(ArQ c 


at 

dx dt 


dt 


dxdt 


dxdt 


(2.300) 


Or A J c = 

/ /(*/■,),.■ 

t=0 X = () 

l f 1 

+ / J (AT 2 )c 

t= 0 x=0 
0 


a / aA 


Ai” + 


a(rAi) a(OA 2 ) 


Af-t 


a.r v a^ 


dt 


at 


+ y^-Y^ix-Xi) 


■=- K. 


ox 2 \ , a(CA 2 ) a(CAi) 


ax \ ax 


i=2 


M-l 


+ 


i=0 
0 r 


(AT,),.- ( 70 — - - 2(7", - y,) 

1 UX 


dt dt 

dt 


+ £ (Ta - W(s - xO 

z=2 


dxdt 


dxdt 


J x=0 


I 

- r 


ox 

(A t 2 )c ( K 2 -^- - 2 (T 2 - y 2 ) 


dt 


J x=0 


£=0 
1 


A A' a(AT,)r 

Aj A] 

Ox 


V 

dt - J 

x-\ ( _ 0 

I 


A 2 A 2 a < A7 ^ 


a.x 


dt 


J X=1 


/ [A 1 C(A7’ 1 )c] t/ dx - j [\ 2 C(AT 2 )c}t f dx 


x=Q 


x=0 
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+ / / Ar 

t Or 0 


(A| — A 2 ) 


d ('/’ 2 - 7 i) 


m 


dx dt 


(2.301) 


Since AJc does not depend (A'i’i)c and (AT 2 )c;, the integrands containing (A'I\)c and 
(A T 2 )c are zero. Thus 


d ( K d\ x \ t d(CAr) 3(CA 2 ) 


(Ac 


(Ac 


+ 


cA 


dt 


d ( K 0X t \ , d(CA 2 ) S(CAi) 


M— 1 

+ £(T, -V,)tf(a: -*<)=() (2.302) 

t=2 

A/ — 1 


03? 

A', 

A 2 


0A, 

dx 
d A 2 


5.x 


+ 


dt dt 

- 2(T\ — Y'i) = 0 at x = 0 

— 2(T 2 — Y 2 ) = 0 at x = 0 


t=2 


dx 

A] — A 2 = 0 at x — 1 
Ai = A 2 = 0 for t = tj 

Equation (2.301) now simplifies to 

'•/ i 


+ E - >' 2 )<5(x - Xi) = o (2.303) 

(2.304) 

(2.305) 

(2.306) 

(2.307) 


A 


Jc = [ f AC 

t=o i=0 


(Ai — A 2 ) 


d(T 2 - TQ 

dt 


dx dt 


From the definition of A J c 


h l 


AJc = J J ACJ'c(x,t.) dx dt 


(=0 1=0 


Comparing Equations (2.308) and (2.309) we have 

, d(T 2 - 1 \ ) 


J'c = (A, - A 2 )- 


dt. 


(2.308) 


(2.309) 


(2.310) 


The inverse technique for reconstructing Ai(2\,2 2 ), A' 2 (7'i/i 2 ) and C{1\ , 7 2 ) can be cast 
in the following algorithmic form: 


1. Assume the properties K x (x,t), A" 2 (x,£) and C(x,i); In the absence of additional 
information, they can be treated as constants. 

2. Generate the corresponding state variables Ti(x, t.) and 2 2 (x, t) by solving Equations 
(2.260-2.267). 
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3. Solve the adjoint problem, Mquations (2.302-2.307), and obtain adjoint variables 
X\(x) and A 2 (. , r) . 

4. Compute' ./'/,• i, ■!' K 2 and .J' r from Equations (2.298), (2.299) and (2.310) respec- 
tively. 

5. Calculate the conjugate coefficients v K \, v K2 and v c form Equations (2.195), (2.190) 
and (2.197) respectively. 

6. Estimate the directions of descent Pki, Pk 2 and Pc from Equations (2.192), (2.193) 
and (2.194) respectively. 

7. Solve the sensitivity problem, Equations (2.268- 2.285) to obtain sensitivity func- 
tions, (A7))/ti, (A7’ 2 )/<i , (A7j)/< 2 , (AT 2 )k 2 , (ATi)( 7 and (A7j)c;. 

8. Compute step sizes ftxi, Pk 2 and /3c using Equations (2.222-2.224 ). 

9. Estimate A'i, K 2 and C using Equations (2.189-2.191 ). 

10. Repeat the above calculation procedure until the discrepancy principle given by 
Equation (2.36) is satisfied. 


2.6 Determination of Constitutive Relationships for 
Oil- Water Flow 


One of the applications of coupled inverse procedure to a practical problem is demon- 
strated in this section. For a flow of two homogeneous immiscible incompressible fluids, 
oil and water for example in an isotropic porous medium, the continuity equations are 
given by [lic.nr( 1 !)72)\ 


dS v , 

0 

7 KI< rw {S w ,S 0 )\ 

dP w 

di 

dx 

A v * ) 

dx 

dS 0 

d 

\(KK r0 {S w ,S 0 )\ 

dP 0 ' 

C di 

dx 

* ) 

' dx 


(2.311) 

(2.312) 
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With the following initial and boundary conditions: 



p,„ 

= Pun 

for 

l = 0 

(2.313) 


Po 

= Poi 

at 

t. = 0 

(2.314) 

KK rw 

dP w 

dx 

— Qw 

at 

x = 0 

(2.315) 

KK rc 

V w 

> dP 0 
dx 

= Qo 

at 

x = 0 

(2.316) 


P W 

= Pwl 

at 

x — L 

(2.317) 


Po 

= Pol 

at 

x = L 

(2.318) 


where the parameters I\ ru , and K r0 are functions of water and oil saturations S w and S 0 . 
Since there are only two fluids in the medium, the sum of S w and S a is unity. Hence K rw 
and K ro are effectively functions of S w . Now we define capillary pressure by 

p — p — p 

1 C — 1 0 1 U) 

Provided near-equilibrium conditions prevail, the capillary pressure, P c is only a function 
of S w . Let water saturation, S w be expressed as 

s w = /(Pc) = f(Po - P w ) 


Etiuafion (2.311) can be written fis 

df(Pc) d 


'KK rw {p c y\ dP w 


Or f 


dt dx 

d£_ dPc_0_ 
d,P r dt dx 


K ) 

dx 


WKI< rw (Pc) 

\ dP w 1 

A 

/ dx 


Or 



(dP 0 

dP w \ 

d 

fKK TW (Pc)) 

dP w 

\ dt 

dt ) 

dx 

A ^ ) 

■ dx 


Or 


WKI< rw (P c )^ 

dP w 

A J 

' dx 


-dP w _ -,0P o 

" dt ° dt 


(2.319) 


where the notation 
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has been used. Since S v , — 1 — S 0l Equation (2.312) can be written as 


OS,, 

Oi 


<)r m 


o_ 

Ox 


KK ro (P c )\ 0P o 


\(KK Tn {P r )\ 

dPo 

LV J 

' Ox 


dx 


c 0Pq _ w 


di 


.dK 

di 


With the following dimensionless quantities 

- P 


:£ l 

L' 

t KPr 

u r L 2 ’ 


w 

T' 


P = K 

0 p: 


c p r ' 


C = e 


d lf_ 

dPr 


K, 


Prw^r 


Kn = 


P ro^r 


we get. dimensionless form of Equations (2.311-2.318) as follows: 

9 ( dP w 
Ox \ w Ox 


± (k ^ 

dx l 0 dx 


-K 
- K 


P w 

Po 

. dP w 


W o 

OX 

op 

° dx 

Pw 


(2.321) 


(2.322) 


C dt 

l 

Q 

CD 2 

(2.323) 

r 0Po 
c dt 

C dt 

(2.324) 

Pwi for 

t = 0 

(2.325) 

Poi for 

t = 0 

(2.326) 

Qw at 

x = 0 

(2.327) 

q a at 

x = 0 

(2.328) 

Pwi at 

x = 1 

(2.329) 

Poi at 

X = 1 

(2.330) 


Equations (2.323-2.330) are similar to Equations (2.260-2.267). Therefore the sensitivity 
Equations and adjoint equations are similar to the sensitivity Equations (2.268-2.285) and 
adjoint Equations (2.302-2.307) respectively derived in Section 2.5. These are repeated 
below for completeness. 


2.6.1 Sensitivity Problem 

Sensitivity Problem for K w 


0 d(AP„) K , 0 ) , 0(, r ,0P„\ „3(A P W ) K ., 0(AP o ) Kw 

& W" Ox ) + & ( AA ” arj ~ c at = “ c at 


(2.331) 
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o ( d(AP f> ) K 
dx \ ° Ox 


n, x _ c d(&P 0 )Kv, = _ c d(A P w ) 


Ot 


in 


d{A P w ) Kw A r , dP w 


-K w v -AK W ^ = 0 at x = 0 

ox ox 


_K 0 ^~ oh<w = 0 at x — 0 

OX 


(AP W ) [< w — {AP 0 ) Kw — 0 9-t •£' — 1 


(AP w ) Kw = (A P 0 )kw = 0 for t = 0 


Sensitivity Problem for K 0 


(2.332) 

(2.333) 

(2.334) 

(2.335) 

(2.336) 


d_ 

dx 


o_ ( u 0(APw) Ko \ Oj APw) Ko _ r 0{AP o ) 

dx v w dx at at 


'k.*&>±) + 1 u K m - c d -^ = - c b{a ^ - k s 

, ox I ox \ ox I ot ot 


d{AP w ) Ko 

-A,„ — — = 0 at x = 0 

OX 


... 0(AP„) K „ . or 0 


■ Ko — -AJC— = 0 at x = Cl 

OX OX 


(A P w ) Ko = (A P 0 ) Ko = 0 at X = 1 


(2.337) 

(2.338) 

(2.339) 

(2.340) 

(2.341) 

(2.342) 


(A P w ) Ko = (A P 0 ) Ko = 0 for t = 0 
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Sensitivity Problem for C 


0 ( ,, i){AP w )c 

OX 


Ox 


K 


dx \ ° Ox 


r ,0{AP w ) c 0P W r ,0{AP o ) 

0 at AC at ~ c at 


(2.343) 

c a(AP„) c , c SPo_ c d[A P,) c 
at. at. at 


(2.344) 

,, 8(A P.)c „ . 

a t« o = o at .r = 0 

ox 


(2.345) 

-K o 0( - Ar ^ C =0 at X = 0 

Ox 


(2.346) 

P„)c = (A P 0 ) c = 0 at x = 1 


(2.347) 

II 

> 

II 

II 


(2.348) 


Adjoint Problem 






K\ ~ + 2 {Pw - P we ) = 0 at X = 0 

OX 


I<2^ + 2 (P 0 - Poe) = 0 at X = 0 


dx 


Aj = A 2 = 0 at x = 1 

A t = A 2 = 0 for t = tj 


(2.351) 

(2.352) 

(2.353) 

(2.354) 
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Gradient. Equations 


•' h w 

J'ko 

J'c 


i)P w i>A, 

dx Ox 
dP 0 d\ 2 



dP 0 

Ot 


(2.355) 

(2.356) 

(2.357) 


2.7 Generalization to Multi-dimensional cases 


The inverse problem discussed in the earlier sections was restricted to one spatial dimen- 
sion. The generalization of this theory is taken up next. The governing equations of 
physical processes can be represented as: 


L{u,pyx,t) = 0 

(2.358) 

where 

u = (ui,u 2 - ■■u N ) r 

(2.359) 

P = {PuP2---Pm) T 

(2.360) 

L = (L u L 2 ---L n ) t 

(2.361) 

x represents the spatial coordinates and t is time. The initial conditions are 


u = fo when t = t 0 

(2.362) 

and the boundary conditions are 

Lub = fi on (Fy/j) 

(2.363) 

Llb = /2 on (Fla) 

(2.364) 


where L UB and L lb are vector operators representing boundary conditions, (u l5 u 2 , ■ ■ • u N ) 
are N state variables and (pi,p 2 , • ■ - Pm) are M parameters. In the direct problem, we need 
to find the unknown state u when parameter p and all initial and boundary conditions are 
known. In contrast to the direct problem, the inverse problem seeks to find the unknown 
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parameters p ( it may also include the unknown initial or boundary conditions) when 
continuous or discrete observations of state u are given. This can be stated as 

find p G Ad 

such that u(p) = z at Ob (2.365) 

where Ad is an admissible set of unknown parameters, and Ob is the observation set of 
space and time where the records of states z are taken. 

Because of measurement and model errors, problem (2.358) has no exact solution. 
The problem 1ms to be solved approximately by an optimization procedure. This is 
equivalent to 


to findp G Ad 

such that J(p) = min J(p) (2.366) 

p€Ad 

where 

t f n jb 

J(p) = I I [ Ui (p) - Zi )] 2 dQ dt (2.367) 

t = o r LB 


where vector [«i(p) — Cj] represents the difference between the model output of u, and its 
observations set Ob. 


2.7.1 Inverse Solution Methods 

Numerical solution of the optimization problem (2.366) can be classified into the following 
three categories: 

• Gauss-Newton 

• Gradient Search 


• Direct Search 
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In order to obtain the Gauss-Newton direction it is necessary to calculate the following 
sensitivity matrix in each iteration of the non-linear least squares minimization: 



thil 

iht'i 

. . . du/ V \ 

Up\ 

dpi 

Op i \ 


du\ 

duz 

<)u N 


OP2 

dpi 

dp 2 


dm 

dU2 

. . . Q u n 

V 

dPM 

dPM 

Op M / 


(2.368) 


where [|^-] is the sensitivity matrix of the discretized state tq with respect to the dis- 
cretized unknown parameter Pj. Gradient search methods are used to avoid the calculation 
of the sensitivity matrix. Hence less computer time is required. Direct search methods do 
not require the calculation of sensitivity matrix (2.368) or the gradient vector but the rate 
of convergence of such methods is generally slow. The grad search method is presented 
in the sections below. 


2.7.2 Direct Problem 

It, is clear that the unknown quantities iq, it 2 • • • U/\r in Equations (2.358-2.364) are assumed 
to be known while px,p 2 , ’"Pm are known. While generating simulated state variables 
u\,u 2 ,- "Un, i.e. pj, p 2 , • • -pM to calculate iq, u 2 , ••• iqv the direct problem (2.358-2.364) 
is nonlinear since the parameters pi,p 2 ," -Pm are functions of tq, u 2 , ■ • ■ u^. Therefore 
an iterative technique is needed for solving the problem a with finite difference method. 

2.7.3 Conjugate Gradient Method For Minimization 

The following iterative process based on the conjugate gradient method is used for esti- 
mation of ]>] , i> 2 ■ ■ - Pm by minimizing the functional J(p). 

pf-'OM) =ft(x,t) - (2.369) 

for n = 0, 1, 2, • • • and i = 1, 2, • • • M where is step sizes for p t in going from iteration 
n to (n + 1) and P? is the direction of descent for p*. The directions of descent are given 

by 




(2.370) 
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which is the conjugation of the gradient direction at the iteration n and the direction 
of descent P" 1 ( .r: , 7. ) at iteration (n — 1). The conjugation coefficient, n" is given by 

1 4 h 


i/' 


I nr;)*dtdx 

I'l f-0 

Y'hjr'ydtdx 

I’l to 


(2.371) 


To perform I, he iteration according to Equation (2.3(h)), we need to compute the step sizes 
/?”,/? 2 , ■ • 'Pm and the gradient of the functional , J ' £, • • • J' n M . In order to determine 
these quantities, the ‘sensitivity problem’ and ‘adjoint problem’ are constructed as follows. 

2.7.4 Sensitivity Analysis 


In order to derive the sensitivity problem for each unknown parameter, we should perturb 
the unknown parameters one at a time. Applying perturbation principle in the Equations 
(2.358,2.362,2.363 and 2.364) the following sensitivity problems are obtained. 


V u L5u - 1 - V p L5p = 0 
5u = 0 when t = t 0 
V u Lc/g<5u -f = 0 on Tub 

Vu-LlbAu + VpL^B^P = 0 on T lb 


(2.372) 

(2.373) 

(2.374) 

(2.375) 


where V„L and V P L are the gradient operators of L with respect to u and p respectively. 
These are given by 
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(2.376) 


(2.377) 


when V u L UIh VpLuB, and V p L lb have the same meaning as above. 

6u is given by the matrix [Auij] where A u i:j denotes the change in u t due to change in 
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the parameter p r The change in parameter pj is given by A pj. Again 8p is given by the 
vector 

8p = [Ap,, Ap 2 ,- ■ ■ Ap M f r 

For the above sensitivity problem we need to change only one parameter at a time. Let 
us consider the effect of change in one parameter p* . In that case all the elements except 
A pk will be zero. For the matrix 8u we have 

A u it j = 0 if j y k 

Auij = A Wjjt if j = k 

The adjoint form of matrices V„L and V P L are given by 
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(2.378) 


(2.379) 


The transpose of the matrices V* and V* are denoted respectively by V„ and V* 
The function to be minimized is given by 

b 

J(p) — J J f[u{p)\dQ,dt (2.380) 


t=to n 


where /[?/,(/;)] is a user-chosen function. [ The choice (T(x)-Yi) 2 is one of the many possible 
for the function /.] The functional J(p) for iteration (n + 1) is obtained as 

b 

J{p n+ ')= j j f[u(p n - (i n P n )]dildl, (2.381) 


t=t 0 n 


Applying Taylor series expansion and neglecting second and higher order terms we get, 

,du 


u{p n - ft n P n ) = u{p n )-(3 n P n 


Op 


u{p n ) - (5 n 8u n 


(2.382) 
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Hence we have' 


•/(/> n+1 ) = f j f\u(p n ) - ft n hv n \d.ildi 

t ’/.» h 


t=t 0 n 


dildt 


(2.383) 


Now, we have to determine f3 n in such a way that J n+1 is minimum. Hence, 

d,J n+l 


dp" 


= 0 


(2.384) 


This will produce M e(|uations for /?", ft%, ■ ■ ■ The (juantities Su n are calculated from 
Equations (2.372-2.375) and P n is estimated from Equations (2.370). From Equations 
(2.370 and 2.371) it is clear that P n is known when the gradient direction J' n is known. 
These gradient directions are calculated from the following adjoint problem. 


2.7.5 Adjoint Problem and Gradient Equations 

We need to find p such that J(p) will be optimized and Equations (2.358-2.364) are 
satisfied. Hence the functional J(p) can be written as, 

tf tf 

J{p)= j J f(p)dttdt+ J J \L{u,p)dttdt (2.385) 

t=t 0 n t=t 0 n 

where A = [A| , A 2 , • • • A^] T are the adjoint functions for the N equations. If N — 1, i.e. for 

one equation model, example heat conduction problem, there is only one adjoint function 

(Ai). 


For a two equations model like oil-water flow through porous media, there are two 
adjoint functions (Ai and A 2 ) 

The variation of the objective function (2.367) 6.J is obtained by perturbing 1 / by 6u in 
Equation (2.367), subtracting from the resulting expression the original Equation (2.367) 
and neglecting the higher order terms. We thus obtain 


h 

" -If 

t=to n 


Of, Of ' 

7T 5u + 5p 
au op 


dildt + j J A [V u L<5u + V p L6p] dildt (2.386) 

t=to n 
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We define the scalar product of any two vector functions u and v on il and the time 
interval [to, tj\ as 


*/ 

(«,u)n iT = J J (uv)diidt (2.387) 

t=t 0 n 

At the boundaries the above product is defined by 

*/ 

( u , v )r UB ,T= [ {uv)\ ruB dt (2.388) 

t=t 0 

b 

( u i v )r LB ,T= j {uv)\ TLB dt (2.389) 

t=t 0 

Using Green’s formula several times to transfer the differential operator from 6u or 5p to 
A we obtain 

f f 

j f X(V n L6u)dttd.t = (A, V u L6u) nt r = (6u, V+LA) n ,r + BIT (2.390) 

b 

j j X{V p L5p)d£ldt'= (A, V p L5p)n,r = ($p, V+LA) n ,r + BIT (2.391) 
«=t 0 n 

where BIT denotes all boundary integral terms and V t |L and V+L are transposed adjoint 
operators of Y„L and V P L respectively. From equations (2.386,2.390 and 2.391 ) we get 
the following expression for 5J: 


SJ 


^<5u, LA + 



+ 


<5p,V±L X + 



+ BIT 
n,T 


where the term BIT in (2.392) is expressed as 


(2.392) 


BIT = (6p,v;L UB X) VuBtT -(5p,W;L LB X) rLB , T + (Su,ViL UB X)r UB ,T 

~ (fob'V+Li jR X)r LB ,T + (<^)iw + {Su^)r LB ,T + {Su,X ) n , t/ - (A ,6u) n ,to 
— {X,XJ u Lij B 5u + V v LuBbp)T UB T ~ u Lib&u + V v Li B bp)r LB ,T (2.393) 


From Equations (2.375-2.375) it is clear that the last three terms of the expression of 
BIT are zero. Again, the variation of the objective function, SJ should not depend 
on the variation of the unknown, Su. Hence the terms containing 5u should be zero. 
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Vanishing the integrands containing 5u we get the following equations for the 

adjoint 

function A. 



V « LX + fu = “ 


(2.394) 

A = 0 

when t = tj 

(2.395) 

o 

II 

■^1 S 3 

+ 

-f- 3 

t> 

on {T un ) 

(2.396) 

-KLlbX+% = 0 

du 

on (14b) 

(2.397) 

Hence Equation (2.392) reduces to 




8J = (fy, VJ LX + |0 ^ + (5p, V+W) rwBt T - (<&P, V+L; jfl A)r LBiT (2.398) 


SI, 7 

Therefore for any unknown parameter p we can obtain the following equation 

'/ 


a/ 

dp 


/ / ( v; 


^ LX 4 — — — 1 dll dt 


t=t 0 n 


dp 


(2.399) 


From the definition of 5J 


8J_ 

dp 


*/ 


— J J J' p dLl dt 


(2.400) 


t=t 0 n 


Comparing Equations (2.399) and (2.400) we have 

■I'v = v;lx + d d (2.401) 

where the boundary integral terms in (2.393) are assumed to be zero. Equation (2.401) is 
called the gradient of the functional with respect to the parameter p. Equations (2.394- 
2.397) is the applicable system for determining the function A. 


2.7.6 Adjoint Operations 

It is possible to find the transpose of the adjoint matrices V*L,Vp L,V„ Lb and V+ Lb by 
following certain symbolic rules. For any matrix A we define A + = (A*) T , where asterisk 
denotes the adjoint operation and superscript T denotes transpose operator. Adjoint 
operation rules for differential operators can be derived from Green’s formula. Some of 
these rules of operation are given in the table below. 
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Table 2.1: Adjoint Operation Rules for Differential Operators 


Element of VL 

Element of V‘L 

F. 

F. 

F Wr 

-U F -) 

&(*) 

F -Wr 

V(FV.) 

V(FV.) 

V(A'V.) 

V(A'V.) 

R(V.) 

-V(V.) 

V(R.) 

-V(V.) 

V(FVG.) 

-FVGV. 

V(ATVG.) 

-AVGV. 

V[AV(F.)] 

FV(A'V.) 


Some of the adjoint operator rules for elements of the boundary operator V Lp are 
given in Table 2.7.6. 

The steps of deriving the adjoint equations can be summarized as follows: 

1. Find gradient operator matrices V u L,V p L,V u Li,B,V u L{/B,VpLpp and V p Lub using 
(2.376) and (2.377). 

2. From (2.378) and (2.379) calculate V*L,V*L,V*Lp and V p *Lp using adjoint oper- 
ation rules, and then obtain V„L,Vp A,V*Lp,V p +Lp by transposition. 

3. Obtain the adjoint Equation (2.394) 

4. Obtain the final and boundary conditions of the adjoint problem (2.39o-2.397). 
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Table 2.2: Adjoint Operation Rules for Elements of the Boundary Operator 


Element, of* V Lb 

Element of V*L W 

(V.)n 

0 

(FV.)n 

(FV.)n 

(A'V.)n 

JS 

E> 

[. KV{F.)]n 

{FI<\ 7> 


5. Obtain all partial derivatives of the cost function J with respect to the parameters 
(2.401). 

2.7.7 Stopping Criterion 

If the problem involves no measurement errors, the traditional condition specified as 

J n+1 {Pi,P2 ■■■Pm) <z (2.402) 

where e is a small specified number, can be used as the stopping criterion. However, the 
recorded data contains measurement error; as a result, the inverse solution will respond to 
the perturbed input data. The estimated parameters consequently will exhibit oscillatory 
behavior as the number of iterations is increased. The computational experience has 
shown that it is advisable to use the discrepancy principle for terminating the iteration 
process in the conjugate gradient method. It is assumed that the standard deviations of 
measurement errors for unknown variables U\,U 2 - are equal in magnitude and. are 
given by a. Then the discrepancy principle that establishes the stopping criterion e can 
be obtained from Equation (2. 307) as 

l ! r ub 

J J a 2 dQ. dt = e 2 

t= 0 I 'm 


(2.403) 
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Then the stopping criterion is taken as 



where e is determined from Equation (2.403). 


(2.404) 


2.7.8 Computational Procedure 

The iterative computational procedure for the conjugate gradient method for the general 
problem can be summarized as follows (assume all the parameters p u p 2 ■ ■ ■ p M are known 
at nth iteration). 

Step 1: Solve the direct problem (2.358-2.364) and compute unknown variables tq, 

«2/ ' ‘ U N- 

Step 2: Solve the adjoint problem (2.394-2.397) and obtain adjoint variables A x , 
A2/ • • A/v. 

Step 3: Compute gradients J'^y, J'* 2 • • -and J' n pM using equation (2.401). 

Step 4: Compute ••• and v™ M from Equation (2.371). Then estimate the 

directions of descent, Pp\,Pp 2 • • • and P" M using Equation (2.370). 

Step 5: Solve the sensitivity problems (2.372-2.375) to obtain 6u. 

Step 6: Compute search step sizes 0p\, (3p 2 • • • and 

Step 7: Compute new set of parameters Pi" +1 ,P 2 n+1 , • • • and p.M n+1 from Equation 
(2.369). 

Step 8: Repeat the above calculational procedure until the discrepancy principle 
defined by Equation (2.404) is satisfied. 



Chapter 3 


SOLUTION OF HEAT 
CONDUCTION EQUATION 


Numerical results obtained using the inverse technique of Chapter 2 are represented below 
for one dimensional heat conduction equations. 


3.1 Opening Remarks 


The accuracy of parameter estimation by the inverse technique depends on variety of 
factors in addition to the measured temperature itself, other quantities of importance are 
time and position measurement, specimen geometry and boundary conditions particularly 
of the heat flux and heat loss variety. In the present study, uncertainties in the measured 
temperature alone have been considered. Other factors that influence the prediction 
include temperature measurement at discrete and data lugging with the finite sampling 
scheme. The fact that the series of direct problems appearing in the inverse procedure 
are also to be numerically solved, contributes to the errors in parameter estimation. 

It has now been established that inverse techniques are mathematically illposed in 
the sense that they amplify errors in the input data. This extreme sensitivity to errors 
manifest in the form of higher errors in predictions as well as possible delay in convergence 
when supplied with an inappropriate initial guess. An unexpected consequence is loss of 
size convergence, namely the absence of convergence in the limit of small space and time 
steps. 
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Ihe results discussed in the present and later chapters have to be interpreted in the 
light of the comments made above. 


3.2 Numerical Issues 


Direct problems arising from the inverse calculation have been solved by a control volume 
finite difference method (Appendix A). The finite difference equations have been solved 
by the TDM A algorithm. When the equation is nonlinear, Picard iterations are required. 

3.3 Estimation of ‘A"’ from Steady State Heat Con- 
duction Problem 

The mathematical formulation for estimating thermal conductivity from steady state heat 
conduction equation has been in Section 2.2. In parameter estimation a detailed examina- 
tion of the sensitivity coefficients can a provide considerable insight into the calculation. 
These coefficients can show possible areas of difficulty and also lead to an improved experi- 
mental design. The sensitivity coefficient is defined as the first derivative of the dependent 
variable with respect to the unknown parameter. If the sensitivity coefficients are either 
small or correlated with one another, the estimation problem becomes very sensitive to 
measurement errors. In Figure (3.1) the distribution of sensitivity of temperature to ther- 
mal conductivity as a function of distance x is shown. The magnitude of the sensitivity 
becomes zero at x = 0 . The sensitivity increases with distance up to x « 0.5. It starts 
decreasing beyond x = 0.5. Finally it goes to zero at x = 1. 

The foregoing results have an implication for the design of the experiment. Since 
the sensitivity is very small at the points close to the boundaries, measurements should 
be taken at the points which are close to the mid-point of the slab. 

To illustrate the validity and accuracy of the conjugate gradient method in predicting 
thermal conductivity K ( T ), we consider a specific example. Here the exact functional form 
of thermal conductivity is assumed to be second-order polynomial with temperature as 
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Figure 3.1: Distribution of sensitivity for K with space for steady state problem 


the dependent, variable, i.e. 

K(T ) = qq ~t~ aixT + a 2 xr 2 (3-1) 

where the constants a 0 , a x and a 2 are taken as 10.0, -0.5 and -0.1 respectively. The two 
boundaries are subjected to constant temperatures, T x = 5.0 and T 2 = 1.0 respectively. 
The grid size is taken as Ax = 0.01 in the finite difference calculations; thermocouple 
spacing Dx equals the finite-difference spacing Ax. 


3.3.1 Inversion of Error-Free Data 

Consider first the inversion of a hypothetical data set assumed to be entirely free of error 
in temperature measurements. The temperature distribution from the solution of the 
direct problem (2.12- 2.14). While simulating measurement temperatures, the form of 
exact K ( T ) is used in the direct problem. Thus the problem is nonlinear and the iterative 
technique is needed for its solution. In the inverse calculation, the thermal conductivity 
takes the form of I< (x ) . The problem thus becomes linear and the estimated temperature 
can be calculated directly. The objective of this section is to show the applicability 
of the inverse procedure in finding K{T ) accurately with no prior information on the 
functional form of the unknown quantities, the so-called function estimation. One of 
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Figure 3.2: Exact and estimated values of K(T) when SD = 0 


the advantages of using the conjugate gradient method is that the initial guesses of the 
unknown quantities can be chosen arbitrarily. This advantage could not be exploited in 
the present study. This is due to the fact that if K{T) satisfies Equation (2.12), aK(T) 
also satisfies Equation (2.12). Here a is a scalar quantity. If the initial guess of K(x) 
is very far from the actual conductivity, the estimated conductivity will differ from the 
actual conductivity. But the estimated and exact functions are parallel to each other. 
In the test case considered here, the initial guess of K(x) used to begin the iteration 
was taken as A'°(x) = 4. The estimated function K (x) obtained using exact input data 
is shown in Figure 3.2. The value of the functional J decreased to less than 10~ 8 after 
iterations, number of iterations is increased. 


3.3.2 Effects of Data Errors 


In order to compare the results for situations involving random measurement errors, we 
assume normally distributed uncorrelated errors with zero mean and a constant standard 
deviation superimposed on the input data. The simulated inexact measurement data Y 
can be expressed as 


Y = Y exact + UJXSD 


(3.2) 
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where 1 'craci fix’ s< >1 1 1 1, i ( >i i of the direct problem with the exact form of I<{T), SD is the 



standard deviation of the measurements, and u is a random variable that is generated 
by subroutine RAND and will be within -2.576 tO 2.576 for a 99% confidence bounds. 
The dimensionless measured temperature with an error corresponding to cr = 0.001 was 
obtained according to Equation (3.2). This represents an error in temperature of about 
0.02%. The inverse solution using these inexact measurements are shown in Figure 3.3. 

Table 3.1: The convergence parameters for steady state heat conduction 


Measurement 
error, a 

Stop 

Criterion 

|j n +! _ jn J 

Number of 
iterations 

CPU 
time (s) 

& rms 

0.000 


10 

0.23 

0.085 


10~ 5 




0.005 


76 

0.52 

0.262 
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3.4 Determination of K and C from a Transient Ex- 
periment 

In the previous section we only the thermal conductivity has been estimated from the 
steady state data. In this section the simultaneous estimation of temperature-dependent 
thermal conductivity and heat capacity has been discussed. No prior information is used 
for the functional forms of the unknown thermal conductivity and heat capacity in the 
present study. 1 wo different boundary conditions have been considered for comparison. 

3.4.1 Non-zero Fluxes at the Boundaries 

The mathematical formulation for simultaneously measuring temperature-dependent ther- 
mal conductivity and heat capacity from a transient heat conduction experiment with 
nonzero fluxes at the boundaries has been presented in Section 2.3. To illustrate the va- 
lidity and accuracy of the conjugate gradient method in simultaneously predicting K(T ) 
and C{T) with inverse analysis from the knowledge of transient temperature records, a 
specific example is considered. Here the exact functional form of thermal conductivity and 
heat capacity are taken as second order polynomials with temperature as the dependent 
variable, i.e. 

K (T) = no + o, x T + u 2 x T 2 (3.3) 

C(T) = b 0 + b l xT + b 2 xT 2 (3.4) 

where the constants a o, a i and a 2 are taken as 6.0, —0.15 and —0.015 respectively, and 
the constants b 0 , b\ and b 2 are 1.2, 0.1 and 0.01 respectively. The medium has initial 
temperature 7’ 0 — 1.0, when t > 0. The two boundaries at x — 0 and x — 1 are subjected 

to constant heat fluxes, q\ = 17.0 and q 2 = 6.0 respectively. To compare the results for 

situations involving random errors, we assume normally distributed uncorrelated errors 
with zero mean and constant standard deviation. The simulated inexact measurement 
data Y(x,t ) can be expressed by 

Y{x,t) = Y exact {x,t)+Lje (3.5) 

While generating simulated measurement temperature V , exact A (1 ) and C(J ) are used 
in the direct problem and thus the problem is nonlinear and the iterative technique in 
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needed for its solution. However, in the inverse calculation, the thermal conductivity and 
heat capacity exist in the form of K(x, t) and C(x,t), so the problem becomes linear and 
the estimated temperature can be calculated directly. 

The space and time increments are taken as Ax = 0.05 and At = 0.02 respectively, 
in the finite difference calculations; the total measurement time is chosen as t f = 1.2. We 
assume that the thermocouples are placed at a regular interval Dx and temperatures are 
measured at these thermocouple points at a regular interval of time Dt. In the entire 
calculations we assume that the thermocouple spacing Dx equals the finite-difference 
spacing Ax. The measurement time step Dt is taken the same as At. 

One of the advantages of using the conjugate gradient method is that the initial 
guesses of the unknown quantities can be chosen arbitrarily. However, this is not valid in 
this problem. The reason is because two unknown functions, K(x,t ) and C(x,t), are to 
be estimated simultaneously by using only the measurement temperature Y(x,t), which 
implies that the estimated temperature T(x,t) obtained by utilizing 'any combination of 
K(x,t) and C(x, /.) could possibly equal Y(x,t ), but the estimated thermal properties are 
not correct ones. 

In order to restrict the region of search directions to obtain the correct inverse 
solutions, a good initial guess of either thermal conductivity, K(x,t), or heat capacity, 
C(x,t), should be given prior to the inverse calculations. Good initial guesses for heat 
capacity can be obtained from from the following energy balance equation 

X—\ 

[q{tj+ 1 ) - q{tj)]At = I C(tj)[T(x, t j+ \) - T{x,tj)]dx (3.6) 

£=0 

where q = q\ — r/ 2 ; j represents the time index; At denotes the time increment for use in 
the finite-difference calculation and C(tj) is an averaged value for heat capacity at t = tj. 

The variation of the sensitivity coefficient of thermal conductivity to temperature as 
a function of time is shown in Figure 3.4 . The sensitivity coefficient is zero at t = 0. The 
sensitivity increases with time, but the sensitivity at x = 1 is smaller compared to that 
at x = 0 for time t ^0. Since the sensitivity of K is small for small time, the estimated 
thermal conductivity is inexact. Therefore it is recommended that the few initial time 
steps be ignored in the inverse calculations. 

The variation of the sensitivity coefficient of heat capacity to temperature as a 
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Figure 3.4: Variation of sensitivity coefficient for K to temperature as a function of time 


function of time is shown in Figure 3.4.1 . The sensitivity is initially zero and its magnitude 
increases with time. Figure 3.4.1 also shows that the sensitivity coefficient does not vary 
with x. It is also observed that the magnitude of the sensitivity coefficient of C is greater 
than that of K. One can then expect reconstruction errors in C to be smaller than in K. 

The foregoing results have several implications for the design of experiments. The 
pattern of sensitivity coefficients indicates that the information obtained from the inverse 
calculations are less reliable for the initial time steps. It also indicates the results for K 
are less critical at the points near x = 0. 

In all the test cases considered here, the initial guesses of K (x, t ) used to begin the 
iteration are taken as 4.5. The estimated functions K(x,t) and C(x,t) with a = 0 V 0, 
o = 0.005 and a = 0.01, arc shown in Figures (3.6-3.13) The value of the functional 
J obtained with a = 0 can be reduced to a very small number as the number of it- 
erations is increased. The dimensionless measured temperature with a = 0.005- and 
u = 0.01 are obtained according to Equation (3.5). In order to show K{T ) and C(T ) 
explicitly as functions of temperature, T, the thermal conductivity and heat capacity at 
z = 0.2, 0.4, 0.6 and 0.8 with a = 0,0.005 and 0.01 are shown in Figures (3.6-3.13). An 
increase in the measurement errors causes a drop in the accuracy of the inverse solution. 
For a = 0.01, the noise level in the reconstructed thermal conductivity function is seen to 
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Figure 3.5: Variation of sensitivity coefficient for C to temperature as a function of time 


be quite high. [Remark: The thermal conductivity reconstruction at a = 0.01 was seen 
to be meaningful after filtering.] 

Figures (3.6-3.13) show that the estimated values of conductivity and heat capacity 
deviate from the actual values for the few initial time steps. The reason is that at interior 
points temperature gradients are initially small. From Equations (2.94-2.95) it is clear 
that if the temperature gradient is small, the gradient of the functional, J'(x,t ) becomes 
small. That implies that the estimated values of the parameters are close to the guess 
values at the points of very low temperature gradients. 

It is to be noted that J’(x,tf) always equals zero since \(x,tf) = 0.0. Therefore the 
estimated values of the unknown functions, K(x , t) and C(x, t) for final time are the same 
as the initial guess values K°(x,tf) and C°(x,tf) respectively. There are two methods 
to avoid such a singularity. One is to use modified conjugate gradient method and the 
second is to record data for a little longer time compared to the actual period of interest, 
hi the present study we applied the second method for solving the gradient equations. 
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Figure 3.7: Exact and estimated values of K{T ) at x — 0.4 
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Figure 3.9: Exact and estimated values of K (T) at x = 0.8 
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Figure 3.13: Exact and estimated values of C(T) at x = 0.8 
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3.4.2 Non-zero Flux at One Boundary and Zero Flux at Oth 

A special case of die flux-flux boundary conditions arises when one of the fluxes is ze 
This confirmation is indeed encpounteied in practice and hence has been seperat 
treated. 2.3.1. 



Figure 3.14: Variation of the sensitiviy coefficient for K with to temperature as a functior. 
of time 


The sensitivity coefficient of the parameters K and C have been shown in the Figures 
(3.14 - 3.15). These figures show low sensitivity coefficients in the region of t = 0. The 
sensitivity become zero at t = 0. Therefore, the information given by inverse solution 
about the parameters will be different from the actual values for the few initial time steps. 
To obtain the actual information we need to ignore a few initial time steps. Figure (3.14) 
also shows that the sensitivity of I< decreases with increasing x. Hence the information 
obtained from the points close to the insulated end is less reliable. The comparison 
between Figures (3.14) and (3.15) shows that the sensitivity to K remains quite, low 
compared to C . Hence the estimated values of C arc expected to be more accurate as 
compared to K. 

To illustrate the validity of the obove statements, we consider a specific example 
where the exact functional form of thermal conductivity and heat capacity are assumed 
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values I< equal to 4.0 and measured temperature with a = 0, a = 0.005 and a = 0.01 
are shown in Figures (3.16 -3.23). As before the functional J obtained for a = 0 can be 
decreased to a small number with increasing number of iterations. Figures (3.16-3 23) 
, shmv that the estimated values of K(T) and C\T) deviate from actual values at the 
region of T = 1.0, he. for the few initial time steps. This is again due to the fact that 
temperature gradient becomes small at some interior points for few initial time steps. In 
Equation 2.82 we see that the adjoint function A is zero at t = t f . This indicates that the 
solution does not improve at the final time, t f . 


Table 3.2: The convergence parameter for the inverse heat conduction with non-zero fluxes 


Measurement 

error, a 

Stop 

Criterion 
| J n+1 - J n \ 

Number of 
iterations 

CPU 
time (s) 

&rms 

r 

v - / rms 

0.000 


97 

12.18 

0.2037 

0.1308 

0.005 

10“ 6 

168 

17.37 

0.2687 

0.1463 

0.010 


105 

12.77 

0.4170 

0.2047 


Table 3.3: The convergence parameter for the inverse heat conduction with non-zero flux 
at one boundary and zero flux at other 


Measurement 
error, a 

Stop 

Criterion 

| jn+i _ jn | 

Number of 
iterations 

CPU 
time (s) 

&rms 

Crms 

0.000 


82 

9.62 

0.4567 

0.1178 

0.005 

10" 6 

70 

8.66 

0.5045 

0.1900 

0.010 


133 

13.23 

0.5744 

0.3630 
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Figure 3.17: Exact and estimated values of K(T ) at x — 0.4 






Figure 3.18: Exact and estimated values of I<(T) at x = 0.6 



Figure 3.19: Exact and estimated values of K(T) at x = 0.8 
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Figure 3.21: Exact and estimated values of C(T) at x 
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Figure 3.23: Exact and estimated values of C(T ) at x = 0.8 
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3.4.3 Flux-Temperature Boundary Conditions 


Tltoro are a lew essential differences between the (lux-flux boundary conditions and (lux- 
temperature boundary conditions, mainly arising from the possibility of a steady state in 
the latter. Inverse calculations with flux-temperature are discussed in the present section. 
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Figure 3.24: Variation of the sensitivity coefficient for K with to temperature as a function 
of time 


The variation of sensitivities of K and C are shown in Figures (3.24 -3.25). The sen- 
sitivities of both K and C are identically zero at x = 1 (where temperature is prescribed). 
This implies that the solutions obtained from the inverse solution will deviate from the 
actual values at the points close to x — 1. The magnitudes of the sensitivities will also 
be in error for a few initial time steps. Therefore, the information given by the inverse 
problem is less reliable for the few initial time steps. The comparison between figures 
(3.24-3.25) shows that the sensitivity of parameter C is about one order of magnitude 
higher than that of K, except near t — 0. This indicates that the solution obtained for C 
is more accurate as compared to K. 

To illustrate the validity of the above statements, we consider a specific example 
where the exact functional form of thermal conductivity and heat capacity are assumed 
to be second-order polynomials with temperature as the dependent variable, i.e. 

K\T) = a„ + a, x T + a 2 / T 2 


(3.9) 
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Figure 3.25: Variation of the sensitivity coefficient for C with to temperature as a function 
of time 


C{T) = 6 0 + bi x T + &2 x T 2 (3.10) 

where the constants a 0 , a x and a 2 for K(T) are taken as 6.0, -0.15 and -0.05 respectively, 
and the constants c 0 , c\ and c 2 for C are chosen as 1.2, 0.15 and 0.05 respectively. The 
initial temperature is assumed to be 1.0. When t > 0, one boundary is subjected to a 
constant flux, q\ = 20.0 and the other one is maintained at a constant temperature of 
T, = 3.0. 

The space and time increments are taken as Ax = 0.025 and At = 0.001 respectively, 
in the finite difference calculations. The total measurement time is chosen as tf = 0.12. 
The thermocouple spacing Dx is assumed to be same as the finite difference spacing Ax. 
The measurement time step Dt is taken as At. The present problem is different from the 
previous two in the sense that the present problem reaches a steady state after some time. 
Using steady state data, C can not be determined from inverse calculation. For the flux- 
flux houndary conditions there is no steady state. In the present calculation, the fluxes 
are not known at the two boundaries, we cannot apply the equation (3.6) to estimate the 
initial guess of C . Therefore, we need to guess both the parameters, K and C arbitrarily. 
However, arbitrary guesses cannot be valid in the present case. The reason is because 
two unknown functions K(x, t) and C(x, t), are to be estimated simultaneously by using 
only the measurement temperature Y(x,t). This implies that the estimated temperature 
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T(x,t) obtained hy utilizing any combination of K(x,t) and C(x,t) could possibly equal 
Y(x,t), but the estimated thermal properties may not the correct ones. It is observed 
from the numerical experiments that the initial guesses of both thermal conductivity and 
heat capacity should be reasonably close to their average values, close to the range. In 
all the test cases considered here, the initial guesses of K(x,t) and C(x,t) used to begin 
the iteration are taken as K°(x,t) = 5.5 and C°(x,t) = 2.0 respectively. 

The thermal conductivity and heat capacity at x = 0.25 and 0.5 with measurement 
errors a = 0, 0.001 and 0.005 are presented in Figures (3.26-3.29). The comparison 
between the figures for estimated K{T ) and C[T) shows that the estimated values of 
heat capacity are more accurate as compared to conductivity. Thus the results agree with 
predictions of the sensitivity analysis. 

Table 3.4: The convergence parameters for flux-temperature boundary conditions 


Measurement 
error, a 

Stop 
Criterion 
|J n+1 - J n | 

Number of 
iterations 

CPU 
time (s) 

I^rms 

Crm5 

0.000 


24 

2.73 

0.3524 

0.5517 

0.005 

10~ 6 

21 

2.53 

0.3549 

0.5518 

0.010 


22 

2.60 

0.3970 

0.5547 




T 

Figure 3.2 7: Exact and estimated values of K(T ) at x = 0.5 











Chapter 4 


COUPLED INVERSE PROBLEM 


4.1 Introduction 

In recent years, researchers and engineers working in the area of ground-water model- 
ing are paying more attention to coupled problem. Examples are flow in saturated and 
unsaturated zones, heat and mass transport in aquifers with and without matrix defor- 
mation and geothermal reservoir evaluation. Coupled problems are also encountered in 
oil reservoirs and nuclear waste repositories. 

The basic concepts and methods of solution for the coupled inverse problem are sim- 
ilar to the single variable inverse problem. However, there are several differences between 
them that make the coupled inverse problem more complex and challenging. First, the 
state variables of a coupled inverse problem often have different dimensions and uncer- 
tainties because they are measured by different instruments. Second, there are crossover 
effects between state variables and parameters. For example, relative permeability of oil 
and water depend upon water saturation, and water saturation depends on capillary pres- 
sure which is a function of oil and water pressures. Finally, there are many options in 
designing an experiment, for the identification of unknown parameters. 

The basic concepts and methods of solution for the coupled inverse problem have 
been discussed systematically in Chapter 2. Numerical results obtained with these tech- 
niques is presented below. 
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1.2 


Estimation of Parameters for Coupled E 
without a Source Term 


Rations 


1.2.1 Steady State Problem 


sensitivity to Parameters 


flip methodology for developing expressions tor use in simultaneously dp. . . 

" ^running un- 

{nown parameters Ad(Tj,T 2 ) and AT 2 (Ti,T 2 ) has been described in Secti 0h 

11 2.4.1. The 

wundary conditions used in this application are of flux-temperature tyn„ rr ,, 

• the distri- 
butions of sensitivity coefficients of K\ and A’ 2 over the domain are shot,,. . 

in r igures 

[4.1-4.21. The magnitude of the sensitivities are seen to decrease with x. fu, 

1 ' y °th the sensi- 
tivities become zero at x = 1. Therefore, the inverse solutions deviate froni „ , . 

' actual values 

at the points those are close to x = 1. 



Figure 4.1: Variation of sensitivity coefficient of K\ with dista^ 


Solution of the Inverse Problem 

To illustrate the validity and establish the accuracy of the conjugate gradj et ^ j n 

simultaneously predicting K l {T 1 ,T 2 ) and A' 2 (T 1 ,7 2 ) by the coupled inve^ analysi$ wg 




4 2 Estimation of Parameters for Coupled Equations without a Source Term 


108 



Figure 4.2: Variation of sensitivity coefficient of I< 2 with distance 


consider a specific example. Here the exact functional form of thermal conductivity and 
heat capacity are assumed to be second-order polynomials with (1\ — T 2 ) as the dependent 
variable, i.e. 

Ki{Ti,T 2 ) = a 0 + «! x (Ti - T 2 ) + a 2 x (T, - T 2 ) 2 (4.1) 

K 2 (T u T 2 ) = b 0 + bi x (Tj — T 2 ) + b 2 x (Ti — T 2 ) 2 (4.2) 

where the constants oo, ui, a 2 and 03 for K\ are taken as 2.0, -0.15 and 0.15 respectively, 
and the constants 6 0 , b\ and b 2 for K 2 are chosen as 3.5, 0.15 and -0.1 respectively. One 
boundary is subjected to constant fluxes q\ and q 2 , while T\ and T 2 are kept constant 

at the other boundary. Here q x and q 2 are specified to be 5.0 and 2.0 respectively. The 

boundary values of both 1\ and T 2 are taken as 1.0. 

The space; increment is taken as Ax = 0.05 in the finite-difference calculations. The 
spacing Dx used for experiments is assumed to same as the finite difference spacing Ax. 

The estimated functions I<\{T\ - T 2 ) and K 2 {T l - T 2 ), obtained when by exact 
measurement data, a — 0.0, are shown in Figures (4.3) and (4.4) respectively. The value 
of the functional .7 obtained in such a case could be decreased to a very small number for 
increasing number of iterations. 

In order to compare the results in the presence of random measurement errors, nor- 
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mally distributed uncorrelated errors with zero moan and a constant, standard deviation 
have been assumed. The measurement, errors for both l\ \ and K 2 are assumed to have 
III,, same standard deviation, n. The simulated inexact, measurement data, Y\ and V a, 
can be expressed as 


Y\ = Y\ exact + Wl<T (4.3) 

b'2 = } 2 exact + UJ 2 0 (4.4) 

where Vi exact fU1 d Y 2(lx act aro the solutions of the direct problem with the exact values of 
Ki{T\ -T 2 ) and K 2 (T\ — T 2 ); ui\ and u> 2 are random variables. The values of these random 
variables are —2.576 to 2.576 with 99% confidence interval. The dimensionless measured 
values Y] and ) 2 with error a = 0.001 are obtained according to Equations (4.3) and 

(4.4). 


From the numerical experiments it has been observed that the initial guesses of the 
parameters, K\ and K 2 should be close to their average values. There are two possible 
reasons for the dependence of the inverse solutions with the initial guesses. First, two 
unknown functions Ki and K 2 , are to be estimated simultaneously by using Y\ and Y 2 . 
This requires that the estimated values T\ and T 2 obtained by utilizing any combination 
of K\ and K 2 could be equal to Y\ and Y 2 , but the estimated properties may not be 
correct. Second, Equations (2.131) and (2.132) show that if K\ and I( 2 satisfy the direct 
problem (2.131- 2.132 ), (\K\ and f3K 2 also satisfy the equations for any scalars a and 
fi. To circumvent these difficulties, the initial guesses of both K\(x) and K 2 (x) used to 
begin the iteration are taken as 3.0. 

The inverse solutions using the inexact measurements are shown in Figures (4.5) 
and (4.6). 
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2.2 Transient Problem 

methodology for simultaneously determining unknown parameters K\ (T\, T 2 ), K 2 (Ti,T 2 ) 
C(T\,T 2 ) has been described in Section 2.4.2. 

irameter Sensitivity 



Figure 4.7: Variation of sensitivity coefficient of K x as a function of time 


Figures (4. 7-4. 9) show the variation of the sensitivity coefficients for K i, K 2 and C 
unction of time. The sensitivity coefficients Ski, Sk 2 and Sc are defined as: 


Ski = 
Si< 2 = 
S c = 


(ATQki ~ (AT 2 )ki 

AKi 

(AXi)k2 ~ (A?2)k2 
A K 2 

(ATi)c - (A T 2 ) c 
| AC 


(4.5) 

‘(4.6) 

(4.7) 


res (4. 7-4. 9) show that the sensitivity coefficients for K i, K 2 and C are very small at 
1 and a maximum at x = 0. The figures also indicates that sensitivities for K 2 and C 
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arc maximum at around t = 0.4. On the other hand, the sensitivity for K\ is maximum 
at t = tf. The above results provide very important information about the solutions of 
the inverse problem. The estimated values of K 2 and C are expected to be close to the 
actual values at the points close to x = 0 and t = 0.4. Similarly, the estimated values of 
K\ will be more accurate for a time close to t — tf. From the adjoint Equation (2.295) 
we see that the adjoint functions are zero at the final time step. Therefore, the estimated 
solutions of K 2 and C will deviate from the actual values near final time step. Figures 
(4. 7-4.9) also show that the sensitivity curves of parameter A 2 are flatter compared to 
the other two parameters. 

Solution of the Inverse Problem 

To illustrate the validity of the conjugate gradient method in simultaneously predicting 
K\{T\,T 2 ), K 2 (T\,T 2 ) and C(T x ,T 2 ) we consider an example. Here the exact functional 
form of A'i, K> and C are taken to be second-order polynomials with ( 7\ — T 2 ) as the 
dependent variable, i.e. 


A'i (1 \ , T 2 ) = a 0 + a,x {T x - T 2 ) + a 2 x (T, - T 2 ) 2 (4.8) 

I< 2 (T x , T 2 ) = b 0 + b x x (Ti - T 2 ) + b 2 x (T x - T 2 ) 2 (4.9) 

C{T X , T 2 ) = c 0 + c, x {Tx - T 2 ) + c 2 x {T x - T 2 ) 2 (4.10) 

where the constants a 0 , a x and a 2 for K\(T\,T 2 ) are taken as 5.0, -1.5 and -0.5 respectively, 

the constants 5 0 , hi and b 2 for K 2 (T\,T 2 ) are taken as 4.5, 1.5 and 1.0 respectively, and 
the constants r 0 , C\ and c 2 for C are chosen as 3.5, -1.5 and -1.0 respectively. The initial 
values of both the state variables T x and T 2 are taken as 1.0. The boundary at x = 0 is 
subjected to constant fluxes q x and q 2 , while at x = 1, Tf and T 2 are kept constant at 
Tu and 7 2 ( resjiectively. In all the test cases considered hare, the values of q x and q 2 are 
taken as 5.0, 1.0 respectively. The values of both r I\i and T 2 i are assumed to be 1.0.* 

The space and time increments are taken as A.r = 0.02 and A t — 0.02 respectively 
in the finite difference calculations. The total measurement time is chosen a stj — 1.2. 
The measurement spacing Dx is assumed to be same as finite difference spacing Ax and 
the measurement time step Dt is taken the same as At of the direct problem. 
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From the numerical experiments it has been observed that the inverse solutions of the 
problem depend on the initial guesses. This is because three unknown functions, K\[x,t) ) 
Kz{x, t) and C(x, t), are to be estimated simultaneously by using only the measured values 
Fi (:/:, t) and Y 2 {x, /,). This requires that the estimated values T \(x, t) and T 2 (x, t) obtained 
by utilizing any combination of Ki(x,t), K 2 (x,t ) and C(x,t ) could possibly be equal to 
Y] (x, t) and Y 2 (. x, t.) respectively. But the estimated properties are not necessarily correct. 
That is why the initial guesses cannot be chosen arbitrarily. The initial guesses should be 
taken to be close to average values of the parameters. If the initial guesses are far from 
the the actual values, the inverse solutions deviate from the actual solutions. In all the 
test cases considered here, initial guesses of Ki{x, t), K 2 (x,t) and C(x,t) used to begin 
the iteration are taken as 5.0, 4.2 and 3.2 respectively. 

The estimated functions of Ki(x,t), K 2 (x,t) and C(:r,f), obtained by using the 
exact measurement data, a — 0.0, are shown in Figures (4.10, 4.12 and 4.14) respectively. 

In order to compare the results for situations involving random measurement errors, 
normally distributed uncorrelated errors with zero mean and constant standard deviation 
are superimposed on the measurement data. The simulated inexact measurement data 
Fi and Y 2 are given by Equations (4.3) and (4.4). In this calculation, it is assumed 
that the standard deviations for both the measured values Y\ and Y 2 are the same. The 
dimensionless measured values of Y\ and Y 2 with errors a = 0.001 are obtained according 
to Equation (4.3) and (4.4). The inverse solutions using these inexact measurements are 
shown in Figures (4.11, 4.13 and 4.15). The figures show the estimated values of the 
parameters at x = 0.25. The estimated values of the parameters K i, K 2 and C have 
been determined in the time interval of t = 0.2 and t = 0.4 since the corresponding errors 
were found to be a minimum. The comparison among the figures show that the estimated 
values of K \ are more accurate as compared to K 2 and C. These results also agree with 
the sensitivity analysis. 















119 


4.3 Estimation of Parameters for Coupled Equations with Source Term 


4.3 Estimation of Parameters for Coupled Equations 
with Source Term 


The mathematical formulations of the problem has been stated in Section 2.5. 

Sensitivity Analysis 



Figure 4.16: Variation of Sensitivity Coefficients of K\ as a function time 


Figures (4.16-4.18) show the distribution of sensitivity coefficients for K i, K 2 and 
C at x = 0.0, 0.5 and 1.0. The sensitivity coefficients Ski, Sk 2 and Sc are given by 
Equations (4.6), (4.7) and (4.7) respectively. The figures show that the sensitivity co- 
efficients are the highest at x = 0 as compared to other locations. The sensitivities of 
the parameters decrease with x and finally go to zero at x = 1. It is observed that the 
sensitivity of K\ at any point x is a maximum at tf. The sensitivity of A 2 is maximum 
for a time close to t = tf. The sensitivity C at any point x is observed to be maximum 
at a time between t = 0.2 and 0.4. The above observations indicate that the information 
obtained from the inverse analysis about the parameter K\ is the best for the points close 
to x = 0 and for time close to final time, tf. From the adjoint Equations (2.295) it is clear 
that the adjoint functions, Ai and A 2 are zero at t = tf Therefore, the inverse analysis 
does not give satisfactory results for the final time step. The above sensitivity analysis 
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indicates that the estimated values of K 2 are more accurate at points close to x = 0 and 
a time close to t, — 0.2. The estimated values of C are expected to he satisfactory at the 
points close to x = 0 for time between t — 0.2 and 0.4. 


Solution of the Inverse Problem 

To illustrate the validity of the above statements, we consider a specific example where 
the exact functional form of K\ and K 2 are assumed to be third-order polynomials while 
C is taken as a second-order polynomial with (T) — T 2 ) as the dependent variable. Hence 

K\ (T, , T 2 ) = go + a, x (Ti - T 2 ) + a 2 x (7j - T 2 ) 2 + a 3 x (T, - T 2 ) 3 (4.11) 

K 2 {T U T 2 ) = h 0 + h, x (7\ - T 2 ) + b 2 x (Ti - T 2 f + b 3 x (T, - T 2 ) 3 (4.12) 

C(T u T 2 ) = c 0 + c, x (2\ — T 2 ) + c 2 x (2\ — 7 2 ) 2 (4.13) 

where the constants a 0 , ai, a 2 and a 3 for Ki(T\,T 2 ) are taken as 5.5, -1.0, -0.5 and -0.1 
respectively, the constants b 0 , 6 X , b 2 and 63 for K 2 (T\,T 2 ) are taken as 5.5, 1.5, 1.0 and 0.5 
respectively, and the constants c 0 , C\ and c 2 are chosen as 1.5, -0.75 and -0.2 respectively. 
The initial values of the variables T\ and T 2 are taken as 1.0 and 2.0 respectively. The 
boundary at x = 0 is subjected to constant fluxes = 6.0 and q 2 — 2.0. The state 
variables 7j and 7 2 are kept at the constant values of 7j; and T 2 i respectively, at the other 
boundary at x = 1. Both the values of 7\ ( and T 2 ; are taken as 1.0. 

The space and time increments are taken as A.t = 0.05 and At = 0.02 respectively, 
in the finite difference calculations. The total measurement time is chosen as tf = 1.0. 
The measurement spacing Dx is assumed to be same as the finite difference spacing Ax 
and the measurement time step Dt is taken as the same as finite difference time step 
At. The estimated functions of Ki(x, t), K 2 (x, t ) and C(x,t ) obtained with measurement 
errors a — 0.0 and o = 0.005 are shown in Figures (4.19-4.24). Figures (4.19), (4.21) and 
(4.23) show the estimated values of the parameters K\, K 2 and C respectively at x = 0.25 
for time between 0.15 and 0.5. The estimated values of the parameters at x = 0.5 for the 
same time range are shown in Figures (4.20), (4.22) and (4.24). 
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Figure 4.20: Exact and estimated values of Ki(Ti — T 2 ) at x 
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Table 1.1: The convergence parameters for coupled steady state problem 


Measurement 
error, a 

Slop 

Criterion 

| jn-fl _ jn j 

Number of 
iterations 

CPU 
time (s) 

(■fti) rms 

(-^2 )rms 

0.000 

10~ 5 

20 

0.16 

0.166 

0.090 

0.001 


20 

0.36 

0.165 

0.159 


Table 4.2: The convergence parameters for coupled inverse problem without source term 


Measurement 
error, a 

Stop 

Criterion 

\ Jn +l _ jn | 

Number of 
iterations 

CPU 
time (s) 

l)rms 

(*&2 )rms 

c rs 

0.000 


142 

30.83 

0.267 

0.58 

1.15 


10 7 






0.001 


128 

28.24 

0.25 

0.59 

1.07 


Table 4.3: The convergence parameters for coupled inverse problem with source term 


Measurement 

error, a 

Stop 

Criterion 

| j n +l _ jn j 

Number of 
iterations 

CPU 
time (s) 

(^l)rms 

(^2 )rms 

C r ms 

0.000 


27 

21.53 

0.385 

0.568 

0.251 


10“ 5 






0.001 


28 

21.55 

0.384 

0.555 

0.332 





Chapter 5 


Oil- Water Flow in an Unsaturated 
Porous Medium 

5.1 Introduction 

1 boro are some important property relationships to be reconstructed in reservoir engi- 
neering calculations. Consider oil-water How in a porous formation, for example. These 
properties include relative permeability of oil, relative permeability of water, capillary 
pressure, porosity of the medium, fluid density, specific heat and thermal conductivity 
of the fluid-saturated rock, viscosity of water and petroleum fluids, enthalpy and specific 
heat for water and petroleum fluids. 

It has been observed that the relative permeability of oil and water and capillary 
pressure depend on water and oil saturation. For homogeneous medium, porosity remains 
constant. The porosity of consolidated rocks (encountered in oil reservoirs) depends on 
the packing of the grains, their shape, arrangement and size distribution. It is obvious 
that the porosity of such medium may vary from point to point. Density of petroleum 
liquids depend on temperature as well. The specific heat and thermal conductivity of fluid 
saturated rock also vary with temperature. The viscosity of oil and water also depend on 
temperature. 

All the above constitutive relationships are extremely important during the simu- 
lation of oil recovery. These parameters are generally not measurable from the physical 
point-of-view. The inverse approach provides a way whereby measurements of state are 
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used to determine the unknown parameters by fitting the model output with the mea- 
surements. 

The present chapter addresses the development of the coupled inverse analysis for 
estimating simultaneously the three parameters, relative permeability of oil, relative per- 
meability of water and capillary pressure, as are encountered in oil-water flow in a porous 
medium. The (low is considered to be isothermal. Also the fluids are assumed to be 
incompressible. 

The mathematical formulation of the inverse problem is given in Section 2.6. 


5.2 Results and Discussions 

The mathematical model of oil-water flow in a porous medium leads the problem of 
coupled equations with Source terms as discussed in Section 4.3. Therefore, the analysis 
of the sensitivity coeliieients for K 0 , K w and C are similar to the sensitivity coefficients for 
K \ , K 2 and C respectively. The parameters K 0 , K w and C are given by Equation (2.322). 
The exact functional form of the parameters K 0 and K w are assumed to be third-order 
polynomials of non-dimensional capillary pressure, P c . The parameter C is taken as a 
second-order polynomial with P r as the dependent variable. The exact functional form of 
these parameters are expressed as 


I<o(Pc) 

== a 0 + a\P c + (i 2 Pc + Pc 

(5.1) 

K w (Pc) 

= 60 + biP c + 62^0 ^ 3 P^ 

(5.2) 

C(P c ) 

= Co + C\P c + C2P c 2 

(5.3) 


where the constants o 0 , cr, a 2 and 03 are taken as 3.2, 1.5, —1.5, and -0.45 respectively, 
the constants b 0 , b\, b 2 and 63 are taken as 3.0, —1.25, 1.25 and 0.4 respectively, and the 
constants c 0 , Ci, c 2 are chosen as 2.0, 1.5 and -0.15, respectively. Oil and water pressure is 
assumed to be constant and equal to 1.0 at t = 0. One boundary is subjected to constant 
fluxes q a = 1 and q w = 1 while the other boundary is kept at constant pressures P a = 1.5 
and P w = 1.0. 


The space and time increments are taken as 0.05 and 0.02 respectively in the finite 
difference calculations. The estimated functions K 0 , K w and C, obtained by using mea- 
surement data with <7 = 0.0, ct = 0.005 are shown in Figures (5. 1-5. 6). Hie estimated 
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parameters are shown for two locations, namely x = 0.2 and x = 0.4. In all the cases, the 
initial guesses ol K„, K w and C are taken as 2.5, 3.0 and 1.5 respectively. 
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Figure 5.1: Exact and estimated values of K W (P C ) at x = 0.2 


^ A 

x ^ A a 

x x A A 

X x A ^ a 

i - \ x A 


EXACT 

SD = 0.0 x 
= 0.005 A 


x 

X x A 


x A a 
x A 


x A 
\ X A 


0.05 0.1 0.15 0.2 0.25 0.3 0.35 

Pr 


Figure 5.2: Exact and estimated values of K W (P C ) at x = 0.4 

















Chapter 6 

Conclusions and Scope for Future 
Work 


6.1 Conclusions 


Details of an inverse technique for a one dimensional nonlinear diffusion problem em- 
ploying the conjugate gradient method have been discussed in the present study. The 
conjugate gradient method provides an efficient, rapidly convergent and a systematic ap- 
proach for the solution of the inverse heat conduction and the coupled inverse problem. 
Numerical results show that the conjugate gradient method does not require a prion 
information on the functional form of the parameters to be estimated. Sensitivity anal- 
yses required during the inverse procedure is seen to be a useful tool for designing an 
experiment. 

The following conclusions have been drawn in the present study: 


Heat Conduction 

1. Parameter estimation is exact if K and C are constants and in pure steady state 
problems. 

2. In general, reconstruction is qualitatively acceptable. 


3. Errors increase with the measurement error, a. 
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4. Larger errors are seen near regions of small temperature gradients and the speci- 
fied temperature boundary. Excessive dependence on the steady state data in an 
unsteady formulation is also seen to contribute to large errors. 

Coupled Diffusion Problem 

1. Reconstruction of the parametric functions is qualitatively acceptable. 

2. Errors increase with the measurement error, a. Errors are generally higher than in 
the single equation model. The convergence rates are sensitive to the initial guesses. 
It is best to prescribe average values of the parameters at the start of the calculation. 

3. Larger errors are seen near regions of small gradients of the dependent variables and 
boundaries with Dirichlet boundary conditions. 

4. The inverse procedure is seen to be applicable for determining constitutive relations 
in oil-water flow in a porous medium. 


6,2 Scope for Future Work 

6.2,1 Experiments 


In all cases considered in the study the input data is taken from numerical results. The 
errors are assumed to be random. In practice, the experimental errors are different. 
Therefore, are needs to estimate the parameters with real experimental data. 


6.2.2 Advection-diffusion Problems 


The inverse procedure should be extended to equations of the type 


C 


'ST 

+ u.VT 

dt 


= V(H(T)VT) + Q 
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Appendix A 


A.l Discretization of the Direct Problem 


Consider the following direct problem wherein the method of discretization of a nonlinear 
equation is presented: 


l = C V)~ sr 1 


(A.l) 


Integrating the above equation over a control volume around the i th node and at n 4* 1 
th time level (Figure A.l), we obtain 
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:r.. i L v 




dx = J \c(T) 
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dr(x, i) 

dt 
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Figure A.l: Control Volume Formulation 
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Or 


'A7V+ 1 + A? +1 \ TP+ 1 - TP +1 


; t+i 


'A t n+1 + A t n _V\ 77+ 1 - 77+ 1 


Ax 


\ 2 ; 

'pra+l _ nrm 

= Cr +1 -- - v 4 Arc 
1 A/, 


Ax 


Or -(ftTf,’ + f<T +l )rr+' + 


(A^ 4 , 1 + 2A7 +1 + A' l n + 1 ) + 2 


,cr\Axi 

At, 


rpn-i-l 


- (A'" +1 + AfJi 1 ) 27L + 1 = 2 qH ^ a ') 2 TP (A.2) 


where the n th level is a known, state and the n + 1 th time level is assumed to be 
unknown. 


A.2 Discretization of the Adjoint Problem and the 
Treatment of the Dirac-Delta Function 


Consider the adjoint problem given by the Equation (2.79): 


3 

i=2 


(A.3) 


dx \7* v ’ ~ J dx ) dt 

Integrating the above equation over a control volume around the i th node and at n — 1 
th time level, we obtain 

X 1 X. . 1 X. . 1 

■- r d [„d\\] n ~ x 

I <te+ / I “ir <tt+ 1 \, 2 
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X. 1 


dx \ ^ dx 


dx + J 

X. 1 

l ~i 


dx+ f 2 ]T [T n ~ l - y n_1 ]5(x - Xi) dx = 0 
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71— 1 
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-i 
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dx + 2(T± - y^- 1 ) = 0 


\ 71. 1 \ 71 1 \Tl-l \7l-l 
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Ax + 2(2?- 1 - y"' 1 ) = 0 
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Or 
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\ 71-1 
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(A. 4) 


In this case, n - 1 th time level is assumed to be the unknown state. The reason is that 
the final conditions are known for the adjoint problem and time marching proceeds from 
t = if to t — 0. 


The same control volume approach has been applied to discretise the adjoint equa- 
tions for the coupled problem as well. 
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